Advanced_scikit_learn_course.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "X_pnHDLAVnRk"
      },
      "source": [
        "# Advanced ML with scikit-learn\n",
        "<p>\n",
        "Bern, 2025<br>\n",
        "Prepared by Dr. Mykhailo Vladymyrov and Matteo Boi.\n",
        "</p>\n",
        "\n",
        "This work is licensed under a <a href=\"https://creativecommons.org/licenses/by-nc-sa/4.0/\">Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License</a>.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HTX7FryF4Fo-"
      },
      "source": [
        "This course is designed to guide students in building complete end to end studies involving Machine learning based tools.\n",
        "\n",
        "\n",
        "* Data inspection with non-linear embedding\n",
        "* Probabililistic gaussian clustering\n",
        "* K means - cluster number selection and metric selection\n",
        "* Clustering discrete data\n",
        "* Data Normalization\n",
        "* Data imputation\n",
        "* Pipelines\n",
        "* Data IO\n",
        "* Visualization techniques and interactive visualizations\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uLOV6Q7s4z9n"
      },
      "source": [
        "We focus here not as much on the tools used, as on the way of thinking, and approaching the task at hand."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "R1kB_-SvVnRp"
      },
      "source": [
        "# Imports and utils"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "x8eD2-bTVnRq"
      },
      "source": [
        "If you work on Colab, you can install the packages with `!pip install package_name`:\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "rQSvgoAicFHQ"
      },
      "outputs": [],
      "source": [
        "!pip install dash hdbscan impyute kmodes pyod umap-learn > /dev/null"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ROBbli_vc-ON"
      },
      "source": [
        "If you work on your local machine, please run your [miniforge](https://github.com/conda-forge/miniforge) or [miniconda](https://docs.conda.io/en/latest/miniconda.html) environment and install the packages like this:\n",
        "\n",
        "```bash\n",
        "conda install dash hdbscan kmodes plotly pyod umap-learn\n",
        "```\n",
        "\n",
        "Some packages may not be installable through conda, but you can usually install them with pip:\n",
        "\n",
        "```bash\n",
        "pip install impyute\n",
        "```\n",
        "\n",
        "To export the environment to a file, you can run:\n",
        "\n",
        "```bash\n",
        "conda env export > environment.yml\n",
        "```\n",
        "\n",
        "While to create the environment from the file:\n",
        "\n",
        "```bash\n",
        "conda env create -f environment.yml\n",
        "```"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "QbxCfwDtVnRq"
      },
      "outputs": [],
      "source": [
        "import os\n",
        "import requests\n",
        "import tarfile\n",
        "\n",
        "import pandas as pd\n",
        "import numpy as np\n",
        "\n",
        "from scipy.stats import kurtosis\n",
        "from scipy.stats import skew\n",
        "\n",
        "from sklearn.model_selection import train_test_split\n",
        "import matplotlib.colors as mcolors\n",
        "\n",
        "# dimensionality reduction\n",
        "import umap\n",
        "from sklearn.decomposition import PCA\n",
        "\n",
        "# clustering kmeans, gaussian\n",
        "from sklearn.cluster import KMeans\n",
        "from sklearn.cluster import HDBSCAN\n",
        "from sklearn.mixture import GaussianMixture\n",
        "from kmodes.kmodes import KModes\n",
        "from kmodes.kprototypes import KPrototypes\n",
        "\n",
        "# metrics for clustering\n",
        "from sklearn.metrics import silhouette_samples, silhouette_score\n",
        "from sklearn.metrics import calinski_harabasz_score\n",
        "from sklearn.metrics import davies_bouldin_score\n",
        "\n",
        "# classification report and metrics\n",
        "from sklearn.metrics import classification_report\n",
        "from sklearn.metrics import accuracy_score, confusion_matrix, f1_score, precision_score, recall_score, roc_auc_score, roc_curve\n",
        "# confusion matrix display\n",
        "from sklearn.metrics import ConfusionMatrixDisplay\n",
        "\n",
        "#regression metrics:\n",
        "from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score\n",
        "\n",
        "# preprocessing\n",
        "from sklearn.preprocessing import PowerTransformer, OrdinalEncoder, OneHotEncoder, StandardScaler, MinMaxScaler, RobustScaler, LabelEncoder\n",
        "\n",
        "# imputations\n",
        "from pyod.models.ecod import ECOD\n",
        "\n",
        "from sklearn.impute import SimpleImputer\n",
        "from sklearn.impute import KNNImputer\n",
        "from sklearn.experimental import enable_iterative_imputer\n",
        "from sklearn.impute import IterativeImputer\n",
        "\n",
        "from impyute.imputation.cs import mice\n",
        "from sklearn.base import BaseEstimator, TransformerMixin\n",
        "\n",
        "# pipelines\n",
        "from sklearn.pipeline import Pipeline\n",
        "from sklearn.compose import ColumnTransformer\n",
        "from sklearn.preprocessing import PolynomialFeatures\n",
        "from sklearn.preprocessing import FunctionTransformer\n",
        "\n",
        "# a model\n",
        "from sklearn.ensemble import RandomForestRegressor\n",
        "from sklearn.model_selection import GridSearchCV\n",
        "\n",
        "# dataset\n",
        "from sklearn.datasets import fetch_california_housing\n",
        "\n",
        "# visualizations\n",
        "import matplotlib.pyplot as plt\n",
        "import seaborn as sns\n",
        "import plotly.express as px\n",
        "from dash import Dash, html, dcc, Input, Output, callback\n",
        "\n",
        "# utils\n",
        "from tqdm.auto import tqdm\n",
        "from joblib import dump, load\n",
        "\n",
        "import warnings\n",
        "warnings.filterwarnings(\"ignore\", category=FutureWarning)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GMSEwruL39o_"
      },
      "source": [
        "# Step 0: study your data"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "6ytLa9qZVnRt"
      },
      "source": [
        "The first step in any data analysis is to study the data. You are the domain expert, and you know your data the best. Yet there are always ways to see the data from a different angle, and to get a better understanding of it.\n",
        "\n",
        "Furthermore, the data might be corrupted, or have missing values, or be in a format that is not suitable for the task at hand.\n",
        "Thus, the first step is to inspect the data yourself, and ensure it is apropriate."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "D7LxTr78VnRt"
      },
      "source": [
        "## First inspection"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "G-UoEIqhVnRt"
      },
      "outputs": [],
      "source": [
        "# load house prices dataset\n",
        "\n",
        "data = fetch_california_housing(as_frame=True) # as_frame=True\n",
        "print(data.keys())"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Gfm5WgMcVnRu"
      },
      "outputs": [],
      "source": [
        "print(data.DESCR)\n",
        "print(data.feature_names)\n",
        "print(data.target_names)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "b3tMmqsqVnRu"
      },
      "outputs": [],
      "source": [
        "x, y = np.array(data.data), np.array(data.target)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "RlWqElXZVnRu"
      },
      "outputs": [],
      "source": [
        "print(x.shape, y.shape)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-1E3ikGFw1gq"
      },
      "source": [
        "As reminder:\n",
        "\n",
        "<img src=\"https://github.com/dsl-unibe-ch/Advanced_ML_with_scikit-learn/raw/main/xy.png\" alt=\"fig_feature\" width=\"110%\"/><br>\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "qRXf9u0kVnRu"
      },
      "outputs": [],
      "source": [
        "df = data.frame\n",
        "df.head()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Tsh3cuSbrOL_"
      },
      "source": [
        "We can now check more into details the data types and if the dataset contains any missing value."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "q8lNlbEarGcR"
      },
      "outputs": [],
      "source": [
        "df.info()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HOHWD-i7tHJq"
      },
      "source": [
        "There is no missing values. Good."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "cKeGosa-tS0W"
      },
      "outputs": [],
      "source": [
        "df.hist(figsize=(12, 10), bins=30, edgecolor=\"black\")\n",
        "plt.subplots_adjust(hspace=0.7, wspace=0.4)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oiqwuFvitcmS"
      },
      "source": [
        "We can start by analyzing features whose distributions align with expectations.  \n",
        "\n",
        "The **median income** shows a distribution with a long tail. This indicates that most incomes follow a roughly normal distribution, but there are a few individuals earning significantly higher salaries.  \n",
        "\n",
        "The **average house age** appears to follow a nearly uniform distribution.  \n",
        "\n",
        "The target variable also has a long-tail distribution. Additionally, there is a threshold effect for high-valued houses—any house priced above 5 is capped at a value of 5.  \n",
        "\n",
        "Looking at **average rooms, average bedrooms, average occupancy**, and **population**, these features span a wide range of values. The bins for the largest values are barely noticeable, suggesting the presence of a few extremely high values that might be outliers. This becomes clearer when we examine the summary statistics for these features."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "liAzgSfRtueg"
      },
      "outputs": [],
      "source": [
        "features_of_interest = [\"AveRooms\", \"AveBedrms\", \"AveOccup\", \"Population\"]\n",
        "df[features_of_interest].describe().T"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "G2BVMKUsxOLA"
      },
      "source": [
        "For each of these features, the large gap between the **maximum** and **75th percentile** values highlights the presence of a few extreme values, confirming our initial intuition.  \n",
        "\n",
        "So far, we haven’t explored **longitude** and **latitude**, which provide geographical information. These features could help us identify whether certain locations are associated with high-priced houses. To do this, we could create a scatter plot with *latitude* and *longitude* on the x- and y-axes, while using the circle size and color to represent house values in each district."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "PJtGQ1xwyAr3"
      },
      "outputs": [],
      "source": [
        "sns.scatterplot(\n",
        "    data=df,\n",
        "    x=\"Longitude\",\n",
        "    y=\"Latitude\",\n",
        "    size=\"MedHouseVal\",\n",
        "    hue=\"MedHouseVal\",\n",
        "    palette=\"viridis\",\n",
        "    alpha=0.5,\n",
        ")\n",
        "plt.legend(title=\"MedHouseVal\", bbox_to_anchor=(1.05, 0.95), loc=\"upper left\")\n",
        "_ = plt.title(\"Median house value depending of\\n their spatial location\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0VL9Sy0myss5"
      },
      "source": [
        "If you’re unfamiliar with California, it’s worth noting that the data points outline the state’s geographical shape. High-priced houses are predominantly located along the coast, where California’s major cities (San Diego, Los Angeles, San Jose, and San Francisco) are situated.\n",
        "\n",
        "We can also create a *pairplot*, which visualizes pairwise relationships between multiple features in a grid of scatter plots. This allows us to observe how features interact with each other and helps identify potential patterns, correlations, or outliers. For example, a pairplot can reveal how house prices relate to variables like median income, house age, or population. It’s particularly useful for quickly spotting trends and understanding the overall structure of the dataset."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "k4ibj4utVnRv"
      },
      "outputs": [],
      "source": [
        "%%time\n",
        "# Drop the unwanted columns\n",
        "columns_drop = [\"Longitude\", \"Latitude\"]\n",
        "subset = df.drop(columns=columns_drop)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "TY3hzh454r3B"
      },
      "outputs": [],
      "source": [
        "# Quantize the target and keep the midpoint for each interval\n",
        "subset[\"MedHouseVal\"] = pd.qcut(subset[\"MedHouseVal\"], 6, retbins=False)\n",
        "subset[\"MedHouseVal\"] = subset[\"MedHouseVal\"].apply(lambda x: x.mid)\n",
        "\n",
        "_ = sns.pairplot(data=subset, hue=\"MedHouseVal\", palette=\"viridis\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "j_2k9Sbw1CB7"
      },
      "source": [
        "In this case, `pd.qcut` divides MedHouseVal into six quantiles, assigning each interval a midpoint for easier plotting and interpretation.\n",
        "\n",
        "Interpreting a pairplot can be challenging due to the large amount of data, but it still offers some valuable insights. We can clearly identify extreme values (potential outliers) in several features. Additionally, it becomes evident that *median income* is a strong indicator for distinguishing high-valued houses from lower-valued ones."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "4E_emOB45NxP"
      },
      "outputs": [],
      "source": [
        "df.plot(kind=\"scatter\", x=\"MedInc\", y=\"MedHouseVal\")\n",
        "plt.xlabel(\"Median Income\")\n",
        "plt.ylabel(\"Median House Value\")\n",
        "plt.title(\"Median Income vs. Median House Value\")\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fUsnLik75M3R"
      },
      "source": [
        "The plot above reveals a strong linear correlation between **Median Income** and **House Value**. Additionally, the House Value capped at 5 appears again, forming distinct horizontal lines. This could pose issues for modeling and should be addressed through preprocessing.\n",
        "\n",
        "Therefore, when building a predictive model, we can expect *longitude*, *latitude*, and *median income* to be important features for predicting *median house values*.\n",
        "\n",
        "However, in the dense areas there are so many points, that we can't really see the distribution. We can use the Kernel Density Estimation (KDE) plot."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fFBx92av7sM1"
      },
      "source": [
        "### Kernel Density Estimation (KDE)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "tnmGp_4O6mMw"
      },
      "outputs": [],
      "source": [
        "# Subsample the data to 10% to reduce computation time\n",
        "subset = df.sample(frac=0.1, random_state=42).drop(columns=columns_drop)\n",
        "print(f\"Subsampled data shape: {subset.shape}\")\n",
        "\n",
        "# Create a PairGrid\n",
        "g = sns.PairGrid(subset, diag_sharey=False, height=2.5)\n",
        "\n",
        "# KDE plot on the upper triangle\n",
        "g.map_upper(sns.kdeplot, cmap='Blues', n_levels=15, fill=True)\n",
        "\n",
        "# Scatter plot on the lower triangle with a regression line\n",
        "g.map_lower(\n",
        "    sns.regplot,\n",
        "    scatter_kws={'color': 'blue', 'alpha': 0.1, 's': 10},\n",
        "    line_kws={'color': 'red', 'alpha': 0.5}\n",
        ")\n",
        "\n",
        "# Histogram on the diagonal\n",
        "g.map_diag(sns.histplot, kde=True, color='royalblue', alpha=0.7, bins=20)\n",
        "\n",
        "# Improve the aesthetics\n",
        "g.fig.suptitle(\"Pairwise Relationships with KDE, Regression, and Histograms\", y=1.02, fontsize=16)\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hP2urSPk8evN"
      },
      "source": [
        "**KDE Plots:** In the upper triangle of the pair grid, we are plotting the Kernel Density Estimation (KDE) to visualize the distribution of features. The KDE plots help us observe the density and underlying structure of the data, which is particularly useful in identifying clusters or regions of higher data concentration. The color gradient in these plots represents the density levels, with darker shades indicating higher concentrations of data points. We set `n_levels=15` to provide a finer resolution of the density contours.\n",
        "\n",
        "**Scatter Plots with Regression Lines:** In the lower triangle, we display scatter plots with regression lines to highlight the relationships between feature pairs. The scatter plot points are semi-transparent (`alpha=0.1`), allowing us to better see areas with high data concentration without overcrowding. The red regression line (`alpha=0.5`) indicates the general trend or linear relationship between the variables.\n",
        "\n",
        "**Histograms with KDE on Diagonal:** On the diagonal, histograms show the distribution of each feature, with a Kernel Density Estimate (KDE) overlaid. The KDE provides a smoother, continuous estimate of the feature distribution, helping us understand how data points are distributed across the feature’s range. The `bins=20` parameter ensures a reasonable level of granularity.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QJX9HA3p73_C"
      },
      "source": [
        "### Principal Component Analysis (PCA)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_8KydH8G83o_"
      },
      "source": [
        "While KDE plots provide us with insight into the density and relationships in the data, they do not necessarily reveal the underlying structure in terms of how the features contribute to variance in the dataset. PCA, on the other hand, helps us reduce the dimensionality of the dataset by transforming the features into a smaller set of components that explain the maximum variance. This can reveal the most important patterns and relationships that may not be immediately apparent in the raw data.\n",
        "\n",
        "Using PCA, we can project the data onto a new set of axes, called **principal components**, which capture the most significant sources of variation in the data. This allows us to visualize the data in lower-dimensional space and identify potential clusters or trends that may be obscured by the complexity of high-dimensional data."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "y-pgZwpUSppk"
      },
      "outputs": [],
      "source": [
        "# Same number of components as the original data\n",
        "n_components = x.shape[1]\n",
        "\n",
        "# Apply PCA to the data\n",
        "pca = PCA(n_components=n_components)  # Reduce to 2 components for 2D visualization\n",
        "x_pca = pca.fit_transform(x)\n",
        "\n",
        "print(x_pca.shape)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ZDaRMuDnT02E"
      },
      "outputs": [],
      "source": [
        "# Explained variance ratio for each principal component\n",
        "explained_variance = pca.explained_variance_ratio_\n",
        "print(f\"Explained Variance Ratio:\\n {explained_variance}\")\n",
        "\n",
        "# Cumulative explained variance\n",
        "cumulative_explained_variance = np.cumsum(explained_variance)\n",
        "print(f\"\\nCumulative Explained Variance:\\n {cumulative_explained_variance}\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "DgHtc3I0VnRv"
      },
      "outputs": [],
      "source": [
        "# Plot the explained variance ratio\n",
        "plt.plot(range(1, len(explained_variance) + 1), explained_variance, color='royalblue', marker='o', linestyle='--', alpha=0.6, label=\"explained variance\")\n",
        "plt.plot(range(1, len(cumulative_explained_variance) + 1), cumulative_explained_variance, color='red', marker='o', linestyle='--', alpha=0.5, label=\"cumulative variance\")\n",
        "plt.title('Expl. Variance Ratio and Cumulative Expl. Variance', fontsize=16)\n",
        "plt.xlabel('Principal Component')\n",
        "plt.ylabel('Explained Variance Ratio')\n",
        "plt.tight_layout()\n",
        "plt.legend()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "tH6sfsfnVnRv"
      },
      "outputs": [],
      "source": [
        "# Plot the PCA results\n",
        "sns.scatterplot(x=x_pca[:, 0], y=x_pca[:, 1], palette='viridis', s=50, alpha=0.7, hue=y)\n",
        "plt.title('PCA of Subsampled Data (2D)')\n",
        "plt.xlabel('Principal Component 1')\n",
        "plt.ylabel('Principal Component 2')\n",
        "plt.grid(True)\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "P3c6LG4pVnRw"
      },
      "outputs": [],
      "source": [
        "# Let's visualize the distribution density with KDE (again, it takes a long time. We could subsample.)\n",
        "\n",
        "sns.kdeplot(x=x_pca[:, 0], y=x_pca[:, 1], cmap='viridis', n_levels=10, fill=True) # n_levels=5, fill=False,\n",
        "plt.title('PCA of Subsampled Data (2D)')\n",
        "plt.xlabel('Principal Component 1')\n",
        "plt.ylabel('Principal Component 2')\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "mKAohH2OX-Q-"
      },
      "outputs": [],
      "source": [
        "# Let's now visualize the third dimension too.\n",
        "\n",
        "fig = px.scatter_3d(x=x_pca[:, 0], y=x_pca[:, 1], z=x_pca[:, 2],\n",
        "                    title=\"3D Scatter Plot of PCA (First 3 Components)\",\n",
        "                    labels={'PC1': 'Principal Component 1',\n",
        "                            'PC2': 'Principal Component 2',\n",
        "                            'PC3': 'Principal Component 3'},\n",
        "                    color=y,\n",
        "                    color_continuous_scale='Viridis',\n",
        "                    opacity=0.7)\n",
        "\n",
        "fig.update_traces(marker=dict(size=3))\n",
        "fig.update_layout(width=1000, height=800)\n",
        "fig.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Qbb8G0zHVnRw"
      },
      "source": [
        "What do you notice about the data, especially in the last two plots?\n",
        "\n",
        "Do you see anything odd?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "NEsdtFWjafKa"
      },
      "outputs": [],
      "source": [
        "df.describe().T"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fMcMYKyNaunp"
      },
      "source": [
        "**Features with Large Variance Dominate**\n",
        "\n",
        "PCA identifies the directions with the highest variance. If some features have much larger values or variances than others (e.g., population), they will dominate the principal components. This causes the PCA to focus primarily on these high-magnitude features, ignoring those with smaller magnitudes."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "12LbBMOMVnRw"
      },
      "source": [
        "## Standardization"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "5aEvT4-jb31m"
      },
      "outputs": [],
      "source": [
        "x_orig = x.copy()\n",
        "y_orig = y.copy()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "o0u8OsEQVnRw"
      },
      "outputs": [],
      "source": [
        "# Standardize the features before applying PCA\n",
        "scaler_x = StandardScaler()\n",
        "scaler_y = StandardScaler()\n",
        "x = scaler_x.fit_transform(x_orig)\n",
        "y = scaler_y.fit_transform(y_orig.reshape(-1, 1)).reshape(-1)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "AJ9XtKR3VnRw"
      },
      "outputs": [],
      "source": [
        "# Print normalization coefficients\n",
        "print(\"Normalization Coefficients:\\n\")\n",
        "\n",
        "print(\"X-axis Feature:\")\n",
        "print(f\"  Mean: {scaler_x.mean_}\")\n",
        "print(f\"\\n  Variance: {scaler_x.var_}\\n\")\n",
        "\n",
        "print(\"Y-axis Feature:\")\n",
        "print(f\"  Mean: {scaler_y.mean_}\")\n",
        "print(f\"  Variance: {scaler_y.var_}\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "shm7SbdajyZm"
      },
      "outputs": [],
      "source": [
        "print(x.var(axis=0)) # variance per features\n",
        "print(y.var())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "z4zrPlFjjrOI"
      },
      "source": [
        "You can always perform the inverse transformation to check if it is correct"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "yDH6hMWcVnRx"
      },
      "outputs": [],
      "source": [
        "x_inv = scaler_x.inverse_transform(x)\n",
        "y_inv = scaler_y.inverse_transform(y.reshape(-1, 1)).reshape(-1)\n",
        "\n",
        "print(x_inv.mean(axis=0), x_inv.var(axis=0))\n",
        "print(y_inv.mean(), y_inv.var())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Ia5XhMyrjvuQ"
      },
      "source": [
        "And compare it with mean and variance of the original data"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "oEkovc0_VnRx"
      },
      "outputs": [],
      "source": [
        "print(x_orig.mean(axis=0), x_orig.var(axis=0))\n",
        "print(y_orig.mean(), y_orig.var())"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "sZTx0I9LfyKp"
      },
      "outputs": [],
      "source": [
        "# PCA on the standardized data\n",
        "pca = PCA(n_components=8)  # Reduce to 2 components for 2D visualization\n",
        "x_pca_standardized = pca.fit_transform(x)\n",
        "\n",
        "sns.scatterplot(x=x_pca_standardized[:, 0], y=x_pca_standardized[:, 1], palette='viridis', s=50, alpha=0.7, hue=y_orig)\n",
        "plt.title('PCA of Subsampled Data (2D)')\n",
        "plt.xlabel('Principal Component 1')\n",
        "plt.ylabel('Principal Component 2')\n",
        "plt.grid(True)\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Explained variance ratio for each principal component\n",
        "explained_variance = pca.explained_variance_ratio_\n",
        "print(f\"Explained Variance Ratio:\\n {explained_variance}\")\n",
        "\n",
        "# Cumulative explained variance\n",
        "cumulative_explained_variance = np.cumsum(explained_variance)\n",
        "print(f\"\\nCumulative Explained Variance:\\n {cumulative_explained_variance}\")\n",
        "\n",
        "# Plot the explained variance ratio\n",
        "plt.plot(range(1, len(explained_variance) + 1), explained_variance, color='royalblue', marker='o', linestyle='--', alpha=0.6, label=\"explained variance\")\n",
        "plt.plot(range(1, len(cumulative_explained_variance) + 1), cumulative_explained_variance, color='red', marker='o', linestyle='--', alpha=0.5, label=\"cumulative variance\")\n",
        "plt.title('Expl. Variance Ratio and Cumulative Expl. Variance', fontsize=16)\n",
        "plt.xlabel('Principal Component')\n",
        "plt.ylabel('Explained Variance Ratio')\n",
        "plt.tight_layout()\n",
        "plt.legend()\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "w6HY2FT29jT6"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "gtxLcEJBcMYY"
      },
      "outputs": [],
      "source": [
        "scaler = MinMaxScaler()\n",
        "x_normalized = scaler.fit_transform(x_orig)\n",
        "\n",
        "# Perform PCA on the normalized data\n",
        "pca_normalized = PCA(n_components=2)\n",
        "x_pca_normalized = pca_normalized.fit_transform(x_normalized)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "_MB9163Xi9Dy"
      },
      "outputs": [],
      "source": [
        "# Plot the normalized PCA results\n",
        "plt.figure(figsize=(12, 6))\n",
        "\n",
        "# Subplot 1: Standardized PCA\n",
        "plt.subplot(1, 2, 1)\n",
        "sns.kdeplot(x=x_pca_standardized[:, 0], y=x_pca_standardized[:, 1], cmap='viridis', n_levels=15, fill=True)\n",
        "plt.title('PCA with Standardization')\n",
        "plt.xlabel('Principal Component 1')\n",
        "plt.ylabel('Principal Component 2')\n",
        "plt.ylim(-5, 5)\n",
        "plt.xlim(-5, 5)\n",
        "\n",
        "# Subplot 2: Normalized PCA\n",
        "plt.subplot(1, 2, 2)\n",
        "sns.kdeplot(x=x_pca_normalized[:, 0], y=x_pca_normalized[:, 1], cmap='viridis', n_levels=15, fill=True)\n",
        "plt.title('PCA with Normalization')\n",
        "plt.xlabel('Principal Component 1')\n",
        "plt.ylabel('Principal Component 2')\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8Omp_FGeg8gl"
      },
      "source": [
        "**Standardization vs Normalization**\n",
        "- Standardization: Centers the data (mean = 0) and scales it to have unit variance (std = 1). Useful for PCA when the data contains features with different scales but similar distributions.\n",
        "\n",
        "- Normalization: Rescales each feature to a range (usually [0, 1]) or divides each sample by its norm (magnitude), making it useful when feature magnitudes vary widely, or when applying PCA to sparse or skewed data.\n",
        "\n",
        "**What it means for PCA**\n",
        "- Standardized PCA: Principal components capture variance by focusing on features with high variance relative to their scale.\n",
        "\n",
        "- Normalized PCA: Each feature contributes more evenly to the principal components, which can help when the dataset contains features with very different units or distributions."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "bBOGGtdohn3f"
      },
      "outputs": [],
      "source": [
        "pd.DataFrame(x, columns=data.data.columns).describe().T"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "_8FCeaGVh2HR"
      },
      "outputs": [],
      "source": [
        "pd.DataFrame(x_normalized, columns=data.data.columns).describe().T"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mBSZ5srZlJCo"
      },
      "source": [
        "### When Standardization is Better (PCA or Linear Regression)  \n",
        "Dataset with features like house size (in m²), number of rooms, and house age (in years). Features have very different scales.\n",
        "\n",
        "PCA and linear regression rely on variance and assume the data is normally distributed. Standardization ensures that all features have equal importance by centering them around zero and scaling them to unit variance. Without standardization, the feature with large range will dominate the analysis, overshadowing other important features.\n",
        "\n",
        "### When Normalization is Better (KNN or Neural Networks)  \n",
        "Image dataset where pixel intensities range from **0 to 255**. You want to train a neural network to classify the images.  \n",
        "\n",
        "Neural networks are sensitive to the scale of the input data. Normalizing pixel intensities to a fixed range (e.g., [0, 1]) speeds up convergence during training and reduces the chance of exploding gradients.  \n",
        "In K-Nearest Neighbors (KNN), normalization ensures that no feature (e.g., pixel brightness) dominates the distance calculation due to its larger magnitude.\n",
        "\n",
        "### **In Short**  \n",
        "- **Standardization** → PCA, linear regression, logistic regression, SVM  \n",
        "- **Normalization** → Neural networks, KNN, image processing, clustering"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "X3BdfWX0VnRx"
      },
      "source": [
        "## Non-linear embedding (UMAP)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_wWYx9hpmY3H"
      },
      "source": [
        "Uniform Manifold Approximation and Projection (UMAP) is a powerful dimensionality reduction technique designed to preserve the structure of high-dimensional data when projecting it to lower dimensions. It excels at maintaining both local and global relationships, making it particularly useful for visualizing complex datasets in 2D or 3D. Unlike PCA, UMAP is non-linear and can capture non-Euclidean structures in data. It’s widely used for tasks like clustering, visualization, and preprocessing for machine learning models. Key parameters, such as the number of neighbors (`n_neighbors`) and minimum distance (`min_dist`), allow fine-tuning the balance between local detail and global structure."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "9R10E_VsVnRx"
      },
      "outputs": [],
      "source": [
        "%%time\n",
        "# Let's generate 2D and 3D umap embeddings of the data\n",
        "# Generate UMAP embeddings\n",
        "umap_2d = umap.UMAP(n_components=2, n_neighbors=30, min_dist=0.1)\n",
        "umap_3d = umap.UMAP(n_components=3, n_neighbors=30, min_dist=0.1)\n",
        "\n",
        "x_umap_2d = umap_2d.fit_transform(x)\n",
        "x_umap_3d = umap_3d.fit_transform(x)\n",
        "\n",
        "print(x_umap_2d.shape, x_umap_3d.shape)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0UU9e77Wmqce"
      },
      "source": [
        "**Note:** The code above is not the same as doing:  \n",
        "```python\n",
        "umap_3d = umap.UMAP(n_components=3, n_neighbors=30, min_dist=0.1)  \n",
        "umap_2d = umap_3d[:, :2]  \n",
        "```\n",
        "\n",
        "Unlike PCA, UMAP does not guarantee that the first two components of a 3D projection will provide the same result as running UMAP directly with `n_components=2`. Each dimensionality reduction is optimized separately and may capture different structures in the data."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "lWU9FJ0QVnRy"
      },
      "outputs": [],
      "source": [
        "# Plot 2D UMAP representation\n",
        "plt.scatter(x_umap_2d[:, 0], x_umap_2d[:, 1], c=y_orig, cmap='viridis', s=15, alpha=0.8)\n",
        "plt.colorbar(label=\"Median House Value\")\n",
        "plt.title(\"UMAP 2D Projection\")\n",
        "plt.xlabel(\"UMAP Component 1\")\n",
        "plt.ylabel(\"UMAP Component 2\")\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "pFnfUNKOVnRy"
      },
      "outputs": [],
      "source": [
        "%%time\n",
        "\n",
        "sns.kdeplot(x=x_umap_2d[:, 0], y=x_umap_2d[:, 1], cmap='viridis', n_levels=15, fill=True)\n",
        "plt.title(\"UMAP 2D Projection\")\n",
        "plt.xlabel(\"UMAP Component 1\")\n",
        "plt.ylabel(\"UMAP Component 2\")\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "mtUPssuAqOFe"
      },
      "outputs": [],
      "source": [
        "import plotly.graph_objects as go\n",
        "\n",
        "# 3D UMAP representation\n",
        "fig = go.Figure(\n",
        "    data=go.Scatter3d(\n",
        "        x=x_umap_3d[:, 0], y=x_umap_3d[:, 1], z=x_umap_3d[:, 2],\n",
        "        mode='markers',\n",
        "        marker=dict(\n",
        "            size=5,\n",
        "            color=y_orig,  # Color by median house value\n",
        "            colorscale='Viridis',\n",
        "            opacity=0.7,\n",
        "            colorbar=dict(title=\"Median House Value\")\n",
        "        )\n",
        "    )\n",
        ")\n",
        "\n",
        "fig.update_layout(\n",
        "    title=\"UMAP 3D Projection\",\n",
        "    scene=dict(\n",
        "        xaxis_title=\"UMAP Component 1\",\n",
        "        yaxis_title=\"UMAP Component 2\",\n",
        "        zaxis_title=\"UMAP Component 3\"\n",
        "    ),\n",
        "    width=1000, height=800\n",
        ")\n",
        "fig.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qMvZQ0fcregU"
      },
      "source": [
        "#### UMAP n_neighbors parameter\n",
        "The **number of neighbors** in UMAP controls the balance between local and global data structure. Smaller values (e.g. 3-15) focus on preserving local patterns, capturing fine-grained relationships between nearby data points. This is useful for detecting clusters or localized features. Larger values (e.g., 30-50) prioritize global structure, offering a broader view of the data by considering more distant neighbors. While this may blur some local details, it helps reveal overarching trends and relationships in the dataset. The choice of `n_neighbors` should depend on whether you are interested in local or global patterns in the data."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "h4zX828uVnRy"
      },
      "outputs": [],
      "source": [
        "%%time\n",
        "# Let's experiment with it using 3, 10, 30, 50, 100, and 300 neighbors.\n",
        "\n",
        "fig, axes = plt.subplots(2, 3, figsize=(15, 10))\n",
        "\n",
        "for i, n_neighbors in enumerate(tqdm([3, 10, 30, 50, 100, 300])):\n",
        "    ax_idx_x = i//3\n",
        "    ax_idx_y = i%3\n",
        "    umap_2d = umap.UMAP(n_components=2, n_neighbors=n_neighbors, min_dist=0.1)\n",
        "    x_umap_2d = umap_2d.fit_transform(x)\n",
        "    sns.scatterplot(x=x_umap_2d[:, 0], y=x_umap_2d[:, 1], hue=y, ax=axes[ax_idx_x, ax_idx_y])\n",
        "\n",
        "    if i > 0:\n",
        "        axes[ax_idx_x, ax_idx_y].get_legend().remove() # this hides the legend for all but the first plot\n",
        "\n",
        "    axes[ax_idx_x, ax_idx_y].set_title(f'n_neighbors={n_neighbors}')\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9zb7WywAVnR2"
      },
      "source": [
        "## Clustering"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WCmkzEmQVnR2"
      },
      "source": [
        "To further analyze the data, we can apply clustering methods to group similar points based on their proximity in the original and in the UMAP space.\n",
        "\n",
        "As with any method we need a way to evaluate the quality. For clustering we can use the silhouette score.\n",
        "The silhouette score is a measure of how similar an object is to its own cluster (cohesion) compared to other clusters (separation).\n",
        "is is defined as follows:\n",
        "\n",
        "$$s(p) = \\frac{b(p) - a(p)}{max(a(p), b(p))}$$\n",
        "\n",
        "where:\n",
        "* $a(p)$ is the average distance of p to all other points in the *same* cluster\n",
        "* $b(p)$ is the average distance of p to all other points in the *closest* cluster\n",
        "* $s(p)$ is the silhouette score of p\n",
        "* $max(a(p), b(p))$ is the maximum of $a(p)$ and $b(p)$\n",
        "\n",
        "The silhouette can be interpreted as follows:\n",
        "* $s(p)$ is between -1 and 1\n",
        "* $s(p)$ is close to 1 if $a(p) << b(p)$, i.e. the point is well clustered\n",
        "* $s(p)$ is close to 0 if $a(p) \\approx b(p)$, i.e. the point is on the border of the cluster\n",
        "* $s(p)$ is close to -1 if $a(p) >> b(p)$, i.e. the point is in the wrong cluster"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### KMeans\n",
        "\n",
        "The basic idea:"
      ],
      "metadata": {
        "id": "5VQxfVicF91_"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "(from https://dashee87.github.io/data%20science/general/Clustering-with-Scikit-with-GIFs/)\n",
        "\n",
        "<img src=\"https://dashee87.github.io/images/kmeans.gif\" alt=\"k-means in action (x marks the spot of the cluster centroid)\" width=\"50%\"/><br>"
      ],
      "metadata": {
        "id": "660-xNbHu4Hb"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "G_Ozyk7xVnR2",
        "scrolled": true
      },
      "outputs": [],
      "source": [
        "# KMeans Clustering on the Original Data with Evaluation Metrics\n",
        "\n",
        "n_clusters_arr = np.arange(2, 11)\n",
        "silhouette_scores = []\n",
        "calinski_harabasz_scores = []\n",
        "davies_bouldin_scores = []\n",
        "\n",
        "for n_clusters in tqdm(n_clusters_arr):\n",
        "    # Set color palette for clusters\n",
        "    cm = plt.colormaps['tab20b']\n",
        "    colors = cm(np.linspace(0, 1, n_clusters))\n",
        "    plt.rcParams['axes.prop_cycle'] = plt.cycler(color=colors)\n",
        "\n",
        "    # Perform KMeans clustering\n",
        "    kmeans = KMeans(n_clusters=n_clusters, n_init=10, max_iter=500, random_state=42)\n",
        "    cluster_labels = kmeans.fit_predict(x)\n",
        "    centroids = kmeans.cluster_centers_\n",
        "\n",
        "    silhouette_vals = silhouette_samples(x, cluster_labels)\n",
        "    mean_silhouette_score = np.mean(silhouette_vals)\n",
        "\n",
        "    fig, axes = plt.subplots(1, 2, figsize=(15, 5))\n",
        "\n",
        "    # Silhouette Plot\n",
        "    y_lower = 0\n",
        "    spacing = int(len(x) / 20)\n",
        "    for i in range(n_clusters):\n",
        "        cluster_vals = silhouette_vals[cluster_labels == i]\n",
        "        cluster_vals.sort()\n",
        "\n",
        "        y_upper = y_lower + len(cluster_vals)\n",
        "        axes[0].barh(\n",
        "            range(y_lower, y_upper), cluster_vals, height=1.0, linewidth=0, label=f'Cluster {i}'\n",
        "        )\n",
        "        cluster_mean = np.mean(cluster_vals)\n",
        "        axes[0].text(0.1, (y_lower + y_upper) / 2, f'{cluster_mean:.3f}', ha='center', va='center')\n",
        "        y_lower = y_upper + spacing\n",
        "\n",
        "    axes[0].axvline(x=mean_silhouette_score, color='red', linestyle='--')\n",
        "    axes[0].set_xlabel('Silhouette Score')\n",
        "    axes[0].set_ylabel('Samples/Cluster')\n",
        "    axes[0].set_xlim([-0.1, 1])\n",
        "    axes[0].set_title('Silhouette Plot')\n",
        "    axes[0].legend(loc='best')\n",
        "\n",
        "    # UMAP Scatter Plot with Cluster Centroids\n",
        "    axes[1].scatter(x_umap_2d[:, 0], x_umap_2d[:, 1], c=cluster_labels, alpha=0.1, cmap=cm)\n",
        "    centroids_umap = umap_2d.transform(centroids)\n",
        "    axes[1].scatter(centroids_umap[:, 0], centroids_umap[:, 1], c='red', s=100, marker='x', label='Centroids')\n",
        "    axes[1].set_title('UMAP with Cluster Centroids')\n",
        "    axes[1].legend()\n",
        "\n",
        "    plt.suptitle(\n",
        "        f'{n_clusters} Clusters | Mean Silhouette: {mean_silhouette_score:.3f} | '\n",
        "        f'Silhouette Score: {silhouette_score(x, cluster_labels):.3f} | '\n",
        "        f'Calinski-Harabasz: {calinski_harabasz_score(x, cluster_labels):.3f} | '\n",
        "        f'Davies-Bouldin: {davies_bouldin_score(x, cluster_labels):.3f}'\n",
        "    )\n",
        "    plt.tight_layout(rect=[0, 0, 1, 0.95])\n",
        "    plt.show()\n",
        "\n",
        "    # Store evaluation metrics\n",
        "    silhouette_scores.append(silhouette_score(x, cluster_labels))\n",
        "    calinski_harabasz_scores.append(calinski_harabasz_score(x, cluster_labels))\n",
        "    davies_bouldin_scores.append(davies_bouldin_score(x, cluster_labels))\n",
        "\n",
        "# Plot Evaluation Metrics for Different Cluster Numbers\n",
        "fig, axes = plt.subplots(1, 3, figsize=(15, 4))\n",
        "axes[0].plot(n_clusters_arr, silhouette_scores, marker='o')\n",
        "axes[0].set(xlabel='Number of Clusters', ylabel='Score', title='Silhouette Score', ylim=[-0.1, 1])\n",
        "\n",
        "axes[1].plot(n_clusters_arr, calinski_harabasz_scores, marker='o')\n",
        "axes[1].set(xlabel='Number of Clusters', title='Calinski-Harabasz Score')\n",
        "\n",
        "axes[2].plot(n_clusters_arr, davies_bouldin_scores, marker='o')\n",
        "axes[2].set(xlabel='Number of Clusters', title='Davies-Bouldin Score')\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "1ljaSdwBVnR3"
      },
      "source": [
        "Silhouette score is a measure of how similar an object is to its own cluster (cohesion) compared to other clusters (separation), higher is better.\n",
        "\n",
        "Calinski Harabasz score is a measure of how dense the clusters are, higher is better.\n",
        "\n",
        "Davies Bouldin score is a measure of how similar the clusters are, lower is better."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "6SmFQPWwVnR3",
        "scrolled": true
      },
      "outputs": [],
      "source": [
        "# KMeans Clustering on the UMAP embeddings with Evaluation Metrics\n",
        "\n",
        "n_clusters_arr = np.arange(2, 11)\n",
        "silhouette_scores = []\n",
        "calinski_harabasz_scores = []\n",
        "davies_bouldin_scores = []\n",
        "\n",
        "for n_clusters in tqdm(n_clusters_arr):\n",
        "    # Set color palette for clusters\n",
        "    cm = plt.colormaps['tab20b']\n",
        "    colors = cm(np.linspace(0, 1, n_clusters))\n",
        "    plt.rcParams['axes.prop_cycle'] = plt.cycler(color=colors)\n",
        "\n",
        "    # Perform KMeans clustering\n",
        "    kmeans = KMeans(n_clusters=n_clusters, n_init=10, max_iter=500, random_state=42)\n",
        "    cluster_labels = kmeans.fit_predict(x_umap_2d)\n",
        "    centroids = kmeans.cluster_centers_\n",
        "\n",
        "    silhouette_vals = silhouette_samples(x_umap_2d, cluster_labels)\n",
        "    mean_silhouette_score = np.mean(silhouette_vals)\n",
        "\n",
        "    fig, axes = plt.subplots(1, 2, figsize=(15, 5))\n",
        "\n",
        "    # Silhouette Plot\n",
        "    y_lower = 0\n",
        "    spacing = int(len(x_umap_2d) / 20)\n",
        "    for i in range(n_clusters):\n",
        "        cluster_vals = silhouette_vals[cluster_labels == i]\n",
        "        cluster_vals.sort()\n",
        "\n",
        "        y_upper = y_lower + len(cluster_vals)\n",
        "        axes[0].barh(\n",
        "            range(y_lower, y_upper), cluster_vals, height=1.0, linewidth=0, label=f'Cluster {i}'\n",
        "        )\n",
        "        cluster_mean = np.mean(cluster_vals)\n",
        "        axes[0].text(0.1, (y_lower + y_upper) / 2, f'{cluster_mean:.3f}', ha='center', va='center')\n",
        "        y_lower = y_upper + spacing\n",
        "\n",
        "    axes[0].axvline(x=mean_silhouette_score, color='red', linestyle='--')\n",
        "    axes[0].set_xlabel('Silhouette Score')\n",
        "    axes[0].set_ylabel('Samples/Cluster')\n",
        "    axes[0].set_xlim([-0.1, 1])\n",
        "    axes[0].set_title('Silhouette Plot')\n",
        "    axes[0].legend(loc='best')\n",
        "\n",
        "    # UMAP Scatter Plot with Cluster Centroids\n",
        "    axes[1].scatter(x_umap_2d[:, 0], x_umap_2d[:, 1], c=cluster_labels, alpha=0.1, cmap=cm)\n",
        "    axes[1].scatter(centroids[:, 0], centroids[:, 1], c='red', s=100, marker='x', label='Centroids')\n",
        "    axes[1].set_title('UMAP with Cluster Centroids')\n",
        "    axes[1].legend()\n",
        "\n",
        "    plt.suptitle(\n",
        "        f'{n_clusters} Clusters | Mean Silhouette: {mean_silhouette_score:.3f} | '\n",
        "        f'Silhouette Score: {silhouette_score(x, cluster_labels):.3f} | '\n",
        "        f'Calinski-Harabasz: {calinski_harabasz_score(x, cluster_labels):.3f} | '\n",
        "        f'Davies-Bouldin: {davies_bouldin_score(x, cluster_labels):.3f}'\n",
        "    )\n",
        "    plt.tight_layout(rect=[0, 0, 1, 0.95])\n",
        "    plt.show()\n",
        "\n",
        "    # Store evaluation metrics\n",
        "    silhouette_scores.append(silhouette_score(x, cluster_labels))\n",
        "    calinski_harabasz_scores.append(calinski_harabasz_score(x, cluster_labels))\n",
        "    davies_bouldin_scores.append(davies_bouldin_score(x, cluster_labels))\n",
        "\n",
        "# Plot Evaluation Metrics for Different Cluster Numbers\n",
        "fig, axes = plt.subplots(1, 3, figsize=(15, 4))\n",
        "axes[0].plot(n_clusters_arr, silhouette_scores, marker='o')\n",
        "axes[0].set(xlabel='Number of Clusters', ylabel='Score', title='Silhouette Score', ylim=[-0.1, 1])\n",
        "\n",
        "axes[1].plot(n_clusters_arr, calinski_harabasz_scores, marker='o')\n",
        "axes[1].set(xlabel='Number of Clusters', title='Calinski-Harabasz Score')\n",
        "\n",
        "axes[2].plot(n_clusters_arr, davies_bouldin_scores, marker='o')\n",
        "axes[2].set(xlabel='Number of Clusters', title='Davies-Bouldin Score')\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### HDBSCAN\n",
        "\n",
        "HDBSCAN (Hierarchical Density-Based Spatial Clustering of Applications with Noise) is an advanced clustering algorithm that builds on DBSCAN, focusing on identifying clusters of varying densities. Unlike KMeans, it doesn't require specifying the number of clusters beforehand and can detect noise (outliers) in the data. HDBSCAN works particularly well for complex datasets where clusters may have irregular shapes or different densities. It constructs a hierarchy of clusters and uses a stability-based approach to select the most meaningful ones. This makes it a robust option for applications like text embeddings or UMAP-reduced data, where traditional clustering methods might struggle."
      ],
      "metadata": {
        "id": "6NSqaXmPGChT"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "vslnE0N0VnR3"
      },
      "outputs": [],
      "source": [
        "%%time\n",
        "# Perform clustering using HDBSCAN on the UMAP-reduced data\n",
        "hdb = HDBSCAN(min_cluster_size=10, n_jobs=2)\n",
        "classes_hdb = hdb.fit_predict(x_umap_2d)\n",
        "\n",
        "# Determine the number of clusters (unique labels)\n",
        "n_clusters = len(np.unique(classes_hdb))\n",
        "\n",
        "cm = plt.colormaps['tab20b']\n",
        "colors = cm(np.linspace(0, 1, n_clusters))\n",
        "plt.rcParams['axes.prop_cycle'] = plt.cycler(color=colors)\n",
        "\n",
        "# Compute silhouette values for the UMAP-reduced data\n",
        "silhouette_values = silhouette_samples(x_umap_2d, classes_hdb)\n",
        "\n",
        "fig, ax = plt.subplots(1,2,figsize=(15,5))\n",
        "\n",
        "n_points = len(x_umap_2d)\n",
        "cluster_y_min = 0\n",
        "ofs = int(n_points / 20)\n",
        "\n",
        "for cluster_idx, cluster_class in enumerate(np.unique(classes_hdb)):\n",
        "\n",
        "    cluster_silhouette_values = silhouette_values[classes_hdb == cluster_class]\n",
        "    cluster_silhouette_values.sort()\n",
        "\n",
        "    cluster_size = len(cluster_silhouette_values)\n",
        "    cluster_y_max = cluster_y_min + cluster_size\n",
        "\n",
        "    ax[0].barh(range(cluster_y_min, cluster_y_max), cluster_silhouette_values, height=1.0, linewidth=0)\n",
        "    ax[0].text(-0.05, (cluster_y_min + cluster_y_max) / 2, str(cluster_class), ha='center', va='center')\n",
        "    cluster_mean_score = np.mean(cluster_silhouette_values)\n",
        "    ax[0].text(0.1, (cluster_y_min + cluster_y_max) / 2, f'{cluster_mean_score:.3f}', ha='center', va='center')\n",
        "\n",
        "    cluster_y_min = cluster_y_max + ofs\n",
        "\n",
        "mean_score = np.mean(silhouette_values)\n",
        "ax[0].axvline(x=mean_score, color='red', linestyle='--')\n",
        "ax[0].set_title(f'Silhouette score')\n",
        "ax[0].set_xlabel('Silhouette score')\n",
        "ax[0].set_ylabel('Sample/Cluster')\n",
        "ax[0].set_yticks([])\n",
        "ax[0].set_xlim([-0.1, 1])\n",
        "\n",
        "# Plot the 2D scatter plot of UMAP-reduced data with HDBSCAN clusters\n",
        "ax[1].scatter(x_umap_2d[:, 0], x_umap_2d[:, 1], c=classes_hdb, alpha=0.1, cmap='tab20b')\n",
        "ax[1].set_title('Scatter plot of the UMAPed data with HDBSCAN clusters')\n",
        "\n",
        "plt.suptitle(f'Silhouette plot for HDBSCAN clusters, mean score: {mean_score:.3f}')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Gaussian clustering"
      ],
      "metadata": {
        "id": "gmiIIEBmaGhN"
      }
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QY2qjDf7VnR4"
      },
      "source": [
        "Gaussian Mixture Models (GMM) are a probabilistic clustering method that assumes the data is generated from a mixture of several Gaussian distributions. Unlike KMeans, which assigns points to the nearest cluster center, GMM provides a **soft clustering approach**: each point is assigned a probability of belonging to each cluster. This makes it more flexible, especially when the clusters have different shapes, sizes, or orientations. GMM is particularly useful when the underlying data distribution is well-modeled by Gaussians, such as in pattern recognition or image segmentation tasks.\n",
        "Let's make clustering of the umapped data."
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Perform GMM clustering on UMAP-reduced data for a range of cluster sizes (2 to 10)\n",
        "\n",
        "n_clusters_arr = np.arange(2, 11)\n",
        "s_s_arr, ch_s_arr, db_s_arr = [], [], []  # Lists to store silhouette, Calinski-Harabasz, and Davies-Bouldin scores\n",
        "\n",
        "for n_clusters in tqdm(n_clusters_arr):\n",
        "\n",
        "    cm = plt.colormaps['tab20b']\n",
        "    colors = cm(np.linspace(0, 1, n_clusters))\n",
        "    plt.rcParams['axes.prop_cycle'] = plt.cycler(color=colors)\n",
        "\n",
        "    # Fit a Gaussian Mixture Model\n",
        "    gmm = GaussianMixture(n_components=n_clusters, n_init=10, max_iter=500)\n",
        "    gmm.fit(x_umap_2d)\n",
        "    classes = gmm.predict(x_umap_2d)\n",
        "    probabilities = gmm.predict_proba(x_umap_2d)\n",
        "\n",
        "    # Extract GMM means and covariances\n",
        "    means = gmm.means_\n",
        "    covariances = gmm.covariances_\n",
        "\n",
        "    # Compute silhouette values\n",
        "    silhouette_values = silhouette_samples(x_umap_2d, classes)\n",
        "\n",
        "    fig, ax = plt.subplots(1, 2, figsize=(15, 5))\n",
        "\n",
        "    n_points = len(x_umap_2d)\n",
        "    cluster_y_min = 0\n",
        "    ofs = int(n_points / 20)\n",
        "\n",
        "    for cluster_idx, cluster_class in enumerate(np.unique(classes)):\n",
        "        cluster_silhouette_values = silhouette_values[classes == cluster_class]\n",
        "        cluster_silhouette_values.sort()\n",
        "\n",
        "        cluster_size = len(cluster_silhouette_values)\n",
        "        cluster_y_max = cluster_y_min + cluster_size\n",
        "\n",
        "        ax[0].barh(range(cluster_y_min, cluster_y_max), cluster_silhouette_values, height=1.0, linewidth=0)\n",
        "        ax[0].text(-0.05, (cluster_y_min + cluster_y_max) / 2, str(cluster_idx), ha='center', va='center')\n",
        "        cluster_mean_score = np.mean(cluster_silhouette_values)\n",
        "        ax[0].text(0.1, (cluster_y_min + cluster_y_max) / 2, f'{cluster_mean_score:.3f}', ha='center', va='center')\n",
        "\n",
        "        cluster_y_min = cluster_y_max + ofs\n",
        "\n",
        "    mean_score = np.mean(silhouette_values)\n",
        "    ax[0].axvline(x=mean_score, color='red', linestyle='--')\n",
        "    ax[0].set_title('Silhouette Score')\n",
        "    ax[0].set_xlabel('Silhouette Score')\n",
        "    ax[0].set_ylabel('Sample/Cluster')\n",
        "    ax[0].set_yticks([])\n",
        "    ax[0].set_xlim([-0.1, 1])\n",
        "\n",
        "    ax[1].scatter(x_umap_2d[:, 0], x_umap_2d[:, 1], c=classes, alpha=0.1, cmap=cm)\n",
        "    ax[1].scatter(means[:, 0], means[:, 1], c='red', s=100, marker='x')  # Mark the centroids\n",
        "    ax[1].set_title('2D Scatter Plot of UMAP Data with Centroids')\n",
        "\n",
        "    plt.suptitle(f'Silhouette Plot for {n_clusters} Clusters, Mean Score: {mean_score:.3f}')\n",
        "    plt.tight_layout()\n",
        "    plt.show()\n",
        "\n",
        "    s_s = silhouette_score(x, classes)\n",
        "    ch_s = calinski_harabasz_score(x, classes)\n",
        "    db_s = davies_bouldin_score(x, classes)\n",
        "\n",
        "    s_s_arr.append(s_s)\n",
        "    ch_s_arr.append(ch_s)\n",
        "    db_s_arr.append(db_s)\n",
        "\n",
        "    fig, ax = plt.subplots(1, n_clusters, figsize=(4 * n_clusters, 4))\n",
        "    cnorm = mcolors.Normalize(vmin=0, vmax=1)\n",
        "    cmap = 'viridis'\n",
        "\n",
        "    for i in range(n_clusters):\n",
        "        ax[i].scatter(x_umap_2d[:, 0], x_umap_2d[:, 1], c=cnorm(probabilities[:, i]), cmap=cmap, alpha=0.1)\n",
        "        ax[i].set(xlabel='UMAP 1', ylabel='UMAP 2', title=f'Cluster {i}')\n",
        "\n",
        "    sm = plt.cm.ScalarMappable(norm=cnorm, cmap=cmap)\n",
        "    fig.colorbar(sm, ax=ax, orientation='horizontal', label='Probability', fraction=0.05)\n",
        "    #plt.tight_layout()\n",
        "    plt.show()\n",
        "\n",
        "# Plot the evaluation scores (silhouette, Calinski-Harabasz, Davies-Bouldin) across the cluster range\n",
        "fig, ax = plt.subplots(1, 3, figsize=(10, 2.5))\n",
        "ax[0].plot(n_clusters_arr, s_s_arr)\n",
        "ax[0].set(xlabel='Number of Clusters', ylabel='Score', title='Silhouette Score', ylim=[-0.1, 1])\n",
        "ax[1].plot(n_clusters_arr, ch_s_arr)\n",
        "ax[1].set(xlabel='Number of Clusters', title='Calinski-Harabasz Score')\n",
        "ax[2].plot(n_clusters_arr, db_s_arr)\n",
        "ax[2].set(xlabel='Number of Clusters', title='Davies-Bouldin Score')\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "zUCnvUZpdkto"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lcMNB6XsVnR4"
      },
      "source": [
        "### Inspect for each class"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "5YA9uHbGVnR5"
      },
      "outputs": [],
      "source": [
        "xy = np.concatenate((x, y.reshape(-1, 1), classes_hdb.reshape(-1, 1)), axis=1)\n",
        "print(xy.shape)\n",
        "\n",
        "df = pd.DataFrame(xy, columns=data.feature_names + data.target_names + ['class_hdbscan'])\n",
        "df.head()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "MD5R8UDkVnR5"
      },
      "outputs": [],
      "source": [
        "%%time\n",
        "# The PairGrid with KDE plots can be slow due to the large dataset size.\n",
        "# To speed up the process, we will subsample the data to 10% of its original size.\n",
        "df_sub = df.sample(frac=0.1)\n",
        "print(f\"Subsampled dataset shape: {df_sub.shape}\")\n",
        "\n",
        "# Create a PairGrid with 'class_hdbscan' as the hue (target variable)\n",
        "g = sns.PairGrid(df_sub, hue='class_hdbscan')\n",
        "\n",
        "# **Upper Triangle: KDE Plot**\n",
        "# We use a KDE (Kernel Density Estimate) plot on the upper triangle to visualize the distribution\n",
        "# of the data in each pair of features. Set the number of contour levels to 10.\n",
        "g.map_upper(sns.kdeplot, n_levels=10)\n",
        "\n",
        "# **Lower Triangle: Scatter Plot with Regression Line**\n",
        "# In the lower triangle, we overlay:\n",
        "# 1. A scatter plot (blue) with low opacity (5%) and small marker size\n",
        "# 2. A regression line (red) with low opacity to complement the scatter plot\n",
        "g.map_lower(sns.regplot,\n",
        "            scatter_kws={'color': 'blue', 'alpha': 0.1, 's': 15},\n",
        "            line_kws={'color': 'red', 'alpha': 0.3})\n",
        "\n",
        "# **Diagonal: Histogram**\n",
        "# On the diagonal, we plot histograms of each feature, stacked for better visualization\n",
        "g.map_diag(sns.histplot, multiple='stack')\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gZ-Jdgs9VnR5"
      },
      "source": [
        "### Discrete KMeans:\n",
        "\n",
        "It's a variation of the traditional KMeans clustering algorithm that is designed to handle categorical or discrete data. Unlike the standard KMeans, which works well with continuous data by minimizing the squared Euclidean distance between points and centroids, discrete KMeans uses a different approach for calculating distances. It typically employs measures like the Hamming distance or other similarity metrics suited for categorical data. The goal remains to partition the data into distinct clusters, but the algorithm adjusts for the unique properties of discrete variables, making it more effective when dealing with categorical features."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "huxfkMoqVnR5"
      },
      "outputs": [],
      "source": [
        "# random categorical data\n",
        "x_test = np.random.choice(20, (1000, 10))\n",
        "\n",
        "km = KModes(n_clusters=4, init='Huang', n_init=10, verbose=1, n_jobs=-1, max_iter=500)\n",
        "\n",
        "clusters = km.fit_predict(x_test)\n",
        "\n",
        "# Print cluster centroids\n",
        "print(km.cluster_centroids_)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ce_WZrT3VnR5"
      },
      "outputs": [],
      "source": [
        "df.columns"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ZTjbssPUVnR5"
      },
      "outputs": [],
      "source": [
        "# We add a column in the original data such that AveRooms_d = int(AveRooms//0.3) and same for AveBedrms_d, Population_d, and HouseAge_d\n",
        "\n",
        "df['AveRooms_d'] = (df['AveRooms'] // 0.3).astype(int)\n",
        "df['AveBedrms_d'] = (df['AveBedrms'] // 0.3).astype(int)\n",
        "df['Population_d'] = (df['Population'] // 0.3).astype(int)\n",
        "df['HouseAge_d'] = (df['HouseAge'] // 0.3).astype(int)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "AQdEUet6VnR5"
      },
      "outputs": [],
      "source": [
        "# Discrete KMeans clustering of the original data\n",
        "\n",
        "x4_d = df[['AveRooms_d', 'AveBedrms_d', 'Population_d', 'HouseAge_d']].values\n",
        "x4 = df[['AveRooms', 'AveBedrms', 'Population', 'HouseAge']].values\n",
        "\n",
        "umap_2d = umap.UMAP(n_components=2, n_neighbors=30, min_dist=0.1)\n",
        "x4_umap_2d = umap_2d.fit_transform(x4)\n",
        "\n",
        "n_clusters_arr = np.arange(2, 6)\n",
        "costs_arr = []\n",
        "\n",
        "for n_clusters in n_clusters_arr:\n",
        "    # Use a more vibrant and distinct palette\n",
        "    palette = sns.color_palette('hsv', n_clusters)  # or 'Set1', 'tab10', 'Paired'\n",
        "\n",
        "    km = KModes(n_clusters=n_clusters, init='Huang', n_init=10, verbose=1, n_jobs=-1, max_iter=500)\n",
        "    clusters = km.fit_predict(x4_d)\n",
        "    costs_arr.append(km.cost_)\n",
        "\n",
        "    # Plot the UMAP-transformed data with clusters\n",
        "    fig, ax = plt.subplots(1, 1, figsize=(6, 6))\n",
        "    sns.scatterplot(\n",
        "        x=x4_umap_2d[:, 0], y=x4_umap_2d[:, 1], hue=clusters, ax=ax,\n",
        "        s=15, alpha=0.7, palette=palette, legend=\"full\"\n",
        "    )\n",
        "    ax.set_title(f'UMAP Scatter Plot with {n_clusters} Clusters')\n",
        "    ax.legend(title='Cluster', bbox_to_anchor=(1.05, 1), loc='upper left')\n",
        "    plt.show()\n",
        "\n",
        "# Plot the cost of clustering\n",
        "plt.figure(figsize=(8, 5))\n",
        "plt.plot(n_clusters_arr, costs_arr, marker='o', color='dodgerblue')\n",
        "plt.xlabel('Number of Clusters')\n",
        "plt.ylabel('Cost')\n",
        "plt.title('Clustering Cost vs Number of Clusters')\n",
        "plt.grid(True)\n",
        "plt.show()\n",
        "plt.close()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Ecy5_xOxVnR6",
        "scrolled": true
      },
      "outputs": [],
      "source": [
        "# You can combine categorical and continuous data in the clustering\n",
        "\n",
        "x8_d = df[['AveRooms_d', 'AveBedrms_d', 'Population_d', 'HouseAge_d', 'MedInc', 'AveOccup', 'Latitude', 'Longitude']].values\n",
        "cat_cols = [0, 1, 2, 3]\n",
        "x8 = df[['AveRooms', 'AveBedrms', 'Population', 'HouseAge']].values\n",
        "\n",
        "umap_2d = umap.UMAP(n_components=2, n_neighbors=30, min_dist=0.1)\n",
        "x8_umap_2d = umap_2d.fit_transform(x8)\n",
        "\n",
        "\n",
        "n_clusters_arr = np.arange(2, 11)\n",
        "costs_arr = []\n",
        "\n",
        "for n_clusters in n_clusters_arr:\n",
        "    palette = sns.color_palette('hsv', n_clusters)\n",
        "\n",
        "    km = KPrototypes(n_clusters=n_clusters, init='Huang', n_init=10, verbose=1, n_jobs=-1, max_iter=500)\n",
        "\n",
        "    clusters = km.fit_predict(x8_d, categorical=cat_cols)\n",
        "    costs_arr.append(km.cost_)\n",
        "\n",
        "    # plot the umap_2d data with the classes\n",
        "    fig, ax = plt.subplots(1, 1, figsize=(3, 3))\n",
        "\n",
        "    sns.scatterplot(x=x8_umap_2d[:, 0], y=x8_umap_2d[:, 1], hue=clusters, ax=ax, s=10, alpha=0.1, palette=palette)\n",
        "    ax.set_title(f'Scatter plot of the UMAPed data with {n_clusters} clusters')\n",
        "    plt.show()\n",
        "\n",
        "plt.plot(n_clusters_arr, costs_arr)\n",
        "plt.xlabel('Number of clusters')\n",
        "plt.ylabel('Cost')\n",
        "plt.title('Cost of clustering')\n",
        "plt.show()\n",
        "plt.close()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9A4DlZgmVnR6"
      },
      "source": [
        "# Step 1: Data Normalization"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "Joqg9qlcVnR6"
      },
      "source": [
        "We briefly discussed the importance of data normalization in the previous chapter.\n",
        "\n",
        "The main reasons for data normalization are:\n",
        "* to make the data more suitable for the model (e.g. neural networks have better convergence properties if the data normalization is compatible with the parameter initialization)\n",
        "* to deal with outliers, skewed distributions, and non-linear relationships\n",
        "* to prevent the model from being biased towards some features\n",
        "* to prevent the model from learnign wrong relationships between features due to improper data representation"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qDZT-enqVnR6"
      },
      "source": [
        "We are going to look at the most standard data preprocessing/normalization techniques:\n",
        "* StandardScaler\n",
        "* MinMaxScaler\n",
        "* RobustScaler\n",
        "* PowerTransformer\n",
        "\n",
        "* LabelEncoder\n",
        "* OrdinalEncoder\n",
        "* OneHotEncoder\n",
        "\n",
        "Which of these do you already know?"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LqYIWjeHVnR6"
      },
      "source": [
        "## Сontinuous data"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "zrrOUldqVnR7"
      },
      "outputs": [],
      "source": [
        "cols = ['AveRooms', 'Longitude', 'MedHouseVal']\n",
        "xy = np.concatenate((x_orig, y_orig.reshape(-1, 1)), axis=1)\n",
        "df_sel = pd.DataFrame(xy, columns=data.feature_names + data.target_names)\n",
        "df_sel = df_sel[cols]\n",
        "\n",
        "v1, v2, v3 = df_sel.values.T\n",
        "\n",
        "g = sns.PairGrid(df_sel)\n",
        "g.map_upper(sns.kdeplot, n_levels=10)\n",
        "g.map_lower(sns.regplot, scatter_kws={'color': 'blue', 'alpha': 0.1, 's': 15}, line_kws={'color': 'red', 'alpha': 0.3})\n",
        "g.map_diag(sns.histplot, multiple='stack')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "lolceEhUVnR7"
      },
      "outputs": [],
      "source": [
        "# For each variable, we perform normalization and plot the distribution (upper row) and the log distribution (lower row)\n",
        "\n",
        "norms = [StandardScaler, MinMaxScaler, RobustScaler, PowerTransformer]\n",
        "norms_names = [ 'StandardScaler', 'MinMaxScaler', 'RobustScaler', 'PowerTransformer']\n",
        "\n",
        "n_norms = len(norms)\n",
        "for v, ttl in zip([v1, v2, v3], cols):\n",
        "    fig, ax = plt.subplots(2, 1 + n_norms, figsize=(15, 4))\n",
        "\n",
        "    for j in range(2):\n",
        "        ax_ji = ax[j, 0]\n",
        "        sns.histplot(v, ax=ax_ji, alpha=0.3)\n",
        "\n",
        "        moments = [v.mean(), v.var(), skew(v), kurtosis(v)]\n",
        "        moments_str = ', '.join([f'{moment:.3f}' for moment, moment_name in zip(moments, ['mean', 'var', 'skew', 'kurt'])])\n",
        "        ax_ji.set_title(f'original data\\n{moments_str}')\n",
        "\n",
        "        if j==1:\n",
        "            ax_ji.set_yscale('log')\n",
        "\n",
        "        for i, (norm, norm_name) in enumerate(zip(norms, norms_names)):\n",
        "            v_norm = norm().fit_transform(v.reshape(-1, 1)).reshape(-1)\n",
        "\n",
        "            # get all first 4 moments of the distribution\n",
        "            moments = [v_norm.mean(), v_norm.var(), skew(v_norm), kurtosis(v_norm)]\n",
        "            moments_str = ', '.join([f'{moment:.3f}' for moment, moment_name in zip(moments, ['mean', 'var', 'skew', 'kurt'])])\n",
        "\n",
        "            ax_ji = ax[j, 1 + i]\n",
        "            sns.histplot(v_norm, ax=ax_ji, alpha=0.3)\n",
        "            ax_ji.set_title(f'{norm_name}\\n{moments_str}')\n",
        "            if j == 1:\n",
        "                ax_ji.set_yscale('log')\n",
        "\n",
        "    plt.suptitle(ttl)\n",
        "    plt.tight_layout(h_pad=2)\n",
        "    plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ONCiDlGCVnR7"
      },
      "source": [
        "As you see the robust scaler is the most robust to outliers, but it's not the best for the other distributions.\n",
        "The power transformer is the best to make data normally distributed, but it's not always producing the best results."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "n_fkfdBzVnR7"
      },
      "source": [
        "## Categorical data"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ZFkMnf4lVnR7"
      },
      "outputs": [],
      "source": [
        "# make dummy dataframe with categorical data:\n",
        "# color: red, green, blue\n",
        "# size: small, medium, large\n",
        "# shape: round, square, triangle, parallelogram, diamond\n",
        "# fill with random values out of the categories, 100 rows\n",
        "\n",
        "n = 100\n",
        "df_cat = pd.DataFrame({\n",
        "    'color': np.random.choice(['red', 'green', 'blue'], n),\n",
        "    'size': np.random.choice([1, 10, 100], n),\n",
        "    'shape': np.random.choice(['round', 'square', 'triangle', 'parallelogram', 'diamond'], n)\n",
        "})\n",
        "\n",
        "\n",
        "n_norms = len(norms)\n",
        "\n",
        "for col in df_cat.columns:\n",
        "    df_col = df_cat[[col]]\n",
        "\n",
        "    normed_vals = LabelEncoder().fit_transform(df_col[col].values)  # output is integer, single column\n",
        "    df_col[f'{col}_LabelEncoder'] = normed_vals\n",
        "\n",
        "    # ensure that the categories are sorted\n",
        "    categories = sorted(np.unique(df_col[col]))\n",
        "    # output is float(integer), multiple columns at a time\n",
        "    normed_vals = OrdinalEncoder(categories=[categories]).fit_transform(df_col[col].values.reshape(-1, 1))\n",
        "    df_col[f'{col}_OrdinalEncoder'] = normed_vals.reshape(-1)\n",
        "\n",
        "    normed_vals = OneHotEncoder(sparse_output=False, drop='first').fit_transform(df_col[col].values.reshape(-1, 1))\n",
        "    df_col[f'{col}_OneHotEncoder'] = [str(normed_vals_i) for normed_vals_i in normed_vals]\n",
        "\n",
        "\n",
        "    print(f'\\ncol={col}')\n",
        "    print(df_col.head(10))\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BeyXkDRQVnR7"
      },
      "source": [
        "So we can see that the OrdinalEncoder is suitable for the size column, as it's naturally ordered.\n",
        "For advanced models label encoding can only be ised if the categories are used as tokens, otherwise it's better to use one hot encoding\n",
        "\n",
        "The one hot can be always used, but it increases the onumber of colums substantially, so it's not always the best choice.\n",
        "To prevent linearly dependent columns, we can use the `drop` parameter, which drops the first column of each category, so we have `n-1` columns for `n` categories."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "82Tdp-PzVnR7"
      },
      "source": [
        "# Step 2: Data Imputation"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "rVYFR3XnVnR7"
      },
      "source": [
        "Data imputation is the process of replacing missing data with substituted values.\n",
        "You often have to deal with missing data, either because the data is not available, or because it's not applicable, or because it's not recorded."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nZT4EuxnVnR8"
      },
      "source": [
        "## Ouliers"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "KnA9ArNfVnR8"
      },
      "outputs": [],
      "source": [
        "ecod = ECOD(n_jobs=-1, contamination=0.1)\n",
        "ecod.fit(x)\n",
        "df['outlier'] = ecod.predict(x)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "dUfcZ9W1VnR8"
      },
      "outputs": [],
      "source": [
        "df['outlier'].value_counts()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "21sNW19gVnR8"
      },
      "outputs": [],
      "source": [
        "# make a scatter plot of the umap_2d data with the outliers\n",
        "fig, ax = plt.subplots(1, 1, figsize=(7, 7))\n",
        "\n",
        "cm = plt.get_cmap('tab20c', 2)\n",
        "plt.rcParams['axes.prop_cycle'] = plt.cycler(color=cm.colors)\n",
        "sns.scatterplot(x=x_umap_2d[:, 0], y=x_umap_2d[:, 1], hue=df['outlier'], ax=ax, s=10, alpha=0.8);"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "B9K-EcuCVnR8"
      },
      "source": [
        "## Missing data"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CDJQjVabVnR8"
      },
      "source": [
        "There several methods of dealing with missing data. The most important is to understand the nature of the missing data:\n",
        "* is it missing completely at random (MCAR) - independent of the observed and missing data\n",
        "* is it missing at random (MAR) - dependent on the observed data but independent of the missing data\n",
        "* is it missing not at random (MNAR) - dependent on the missing data\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "P6Ejyzr_VnR8"
      },
      "source": [
        " [Summary-of-methods-for-handling-missing-values-in-data-instances](https://www.researchgate.net/figure/Summary-of-methods-for-handling-missing-values-in-data-instances-259_fig12_313510665)\n",
        "\n",
        "\n",
        "<img src=\"https://github.com/dsl-unibe-ch/Advanced_ML_with_scikit-learn/raw/main/fig_imputation.png\" alt=\"fig_imputation\" width=\"50%\"/><br>"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "FSxk3ebKVnR8"
      },
      "outputs": [],
      "source": [
        "# make a copy of the original data\n",
        "df_imp = df.copy()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "0t8A39vlVnR9"
      },
      "outputs": [],
      "source": [
        "xy = np.concatenate((x, y.reshape(-1, 1)), axis=1)\n",
        "print(xy.shape)\n",
        "df_imp = pd.DataFrame(xy, columns=data.feature_names + data.target_names)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "3VUCMS14VnR9"
      },
      "outputs": [],
      "source": [
        "# make some missing values in the data\n",
        "df_imp.loc[df_imp.sample(frac=0.5).index, 'AveRooms'] = np.nan"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "5U6CkDHwVnR9"
      },
      "outputs": [],
      "source": [
        "df_imp.isna().mean()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "poYL0bv3VnR9"
      },
      "outputs": [],
      "source": [
        "# We now impute the missing values with the mean\n",
        "imp = SimpleImputer(strategy='mean')  # strategy can be 'mean', 'median', 'most_frequent', 'constant' (according to `fill_value`)\n",
        "\n",
        "df_imp['AveRooms_imp'] = imp.fit_transform(df_imp[['AveRooms']])\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "-iwSUCgKVnR9"
      },
      "outputs": [],
      "source": [
        "# Distribution plot of the original and imputed data\n",
        "\n",
        "fig, ax = plt.subplots(1, 1, figsize=(15, 4))\n",
        "\n",
        "bins=np.linspace(-10, 30, 40)\n",
        "# alpha is .3 to make the plot more transparent, log scale to see the small values\n",
        "sns.histplot(df['AveRooms'], ax=ax, label='Original', alpha=0.3, color='red', bins=bins)\n",
        "sns.histplot(df_imp['AveRooms_imp'], ax=ax, label='Imputed', alpha=0.3, color='blue', bins=bins)\n",
        "ax.set_yscale('log')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "49DWNoaoVnR9"
      },
      "outputs": [],
      "source": [
        "df_imp.columns"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "uHnYPwUEVnR9"
      },
      "outputs": [],
      "source": [
        "knn = KNNImputer(n_neighbors=5)\n",
        "\n",
        "imp_cols = ['MedInc', 'HouseAge', 'AveRooms', 'AveBedrms', 'Population', 'AveOccup', 'Latitude', 'Longitude']\n",
        "df_imp['AveRooms_imp_knn'] = knn.fit_transform(df_imp[imp_cols])[:, 2]\n",
        "\n",
        "sns.histplot(df_imp['AveRooms_imp_knn'][df_imp['AveRooms'].isna()], label='Imputed KNN', alpha=0.3);"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "mMkcuTGSVnR9"
      },
      "outputs": [],
      "source": [
        "# Distribution plot of the original and imputed data\n",
        "\n",
        "fig, ax = plt.subplots(1, 1, figsize=(15, 4))\n",
        "\n",
        "bins=np.linspace(-10, 30, 40)\n",
        "# alpha is .3 to make the plot more transparent, log scale to see the small values\n",
        "sns.histplot(df['AveRooms'], ax=ax, label='Original', alpha=0.3, color='red', bins=bins)\n",
        "# sns.histplot(df_imp['AveRooms_imp'], ax=ax, label='Imputed', alpha=0.3, color='blue', bins=bins)\n",
        "sns.histplot(df_imp['AveRooms_imp_knn'], ax=ax, label='Imputed KNN', alpha=0.3, color='green', bins=bins)\n",
        "ax.set_yscale('log')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "UOqWGCXrVnR-"
      },
      "outputs": [],
      "source": [
        "# Iterative imputer of the data:\n",
        "# Multivariate imputer that estimates each feature from all the others.\n",
        "imp = IterativeImputer()\n",
        "\n",
        "df_imp['AveRooms_imp_iter'] = imp.fit_transform(df_imp[imp_cols])[:, 2]\n",
        "sns.histplot(df_imp['AveRooms_imp_iter'][df_imp['AveRooms'].isna()], label='Imputed regression', alpha=0.3);"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "_L3w_ZrgVnR-"
      },
      "outputs": [],
      "source": [
        "# Distribution plot of the original and imputed data\n",
        "\n",
        "fig, ax = plt.subplots(1, 1, figsize=(15, 4))\n",
        "\n",
        "bins=np.linspace(-10, 30, 40)\n",
        "# alpha is .3 to make the plot more transparent, log scale to see the small values\n",
        "sns.histplot(df['AveRooms'], ax=ax, label='Original', alpha=0.3, color='red', bins=bins)\n",
        "#sns.histplot(df_imp['AveRooms_imp'], ax=ax, label='Imputed', alpha=0.3, color='blue', bins=bins)\n",
        "#sns.histplot(df_imp['AveRooms_imp_knn'], ax=ax, label='Imputed KNN', alpha=0.3, color='green', bins=bins)\n",
        "sns.histplot(df_imp['AveRooms_imp_iter'], ax=ax, label='Imputed Iter', alpha=0.3, color='yellow', bins=bins)\n",
        "ax.set_yscale('log')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "bSvbHZY4VnR-"
      },
      "outputs": [],
      "source": [
        "# np.float = np.float64  # might be needed for new np for impyute to run"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "df_imp['AveRooms_imp_mice'] = mice(df_imp[imp_cols].values)[:, 2]\n",
        "\n",
        "sns.histplot(df_imp['AveRooms_imp_mice'][df_imp['AveRooms'].isna()], label='Imputed MICE', alpha=0.3);"
      ],
      "metadata": {
        "id": "KVI8ry2nb95t"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "nHuj0zwXVnR-"
      },
      "outputs": [],
      "source": [
        "# Distribution plot of the original and imputed data\n",
        "\n",
        "fig, ax = plt.subplots(1, 1, figsize=(15, 4))\n",
        "\n",
        "bins=np.linspace(-10, 30, 40)\n",
        "# alpha is .3 to make the plot more transparent, log scale to see the small values\n",
        "sns.histplot(df['AveRooms'], ax=ax, label='Original', alpha=0.3, color='red', bins=bins)\n",
        "#sns.histplot(df_imp['AveRooms_imp'], ax=ax, label='Imputed', alpha=0.3, color='blue', bins=bins)\n",
        "#sns.histplot(df_imp['AveRooms_imp_knn'], ax=ax, label='Imputed KNN', alpha=0.3, color='green', bins=bins)\n",
        "sns.histplot(df_imp['AveRooms_imp_iter'], ax=ax, label='Imputed Iter', alpha=0.3, color='yellow', bins=bins)\n",
        "sns.histplot(df_imp['AveRooms_imp_mice'][df_imp['AveRooms'].isna()], label='Imputed MICE', alpha=0.3, color='purple', bins=bins)\n",
        "ax.set_yscale('log')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GJM3J9pVVnR-"
      },
      "source": [
        "Please explore the other imputation methods in the `impyute` library.\n",
        "There is a number of them for different types of data."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wJ5Duc1IVnR-"
      },
      "source": [
        "# Step 3. Data Transformation Pipeline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "cWZHklo4VnR_"
      },
      "source": [
        "We see by now that we have a number of preprocessing steps even before we start fitting the model. To simplify our work we can use the `Pipeline` class from `sklearn.pipeline` module.\n",
        "\n",
        "The aim of a pipeline is to assemble several preprocessing steps and a model in a single object."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "X9q-khK_VnR_"
      },
      "outputs": [],
      "source": [
        "data = fetch_california_housing()\n",
        "x, y = data.data, data.target\n",
        "\n",
        "col_names = data.feature_names\n",
        "tgt_names = data.target_names\n",
        "\n",
        "xy = np.concatenate((x, y.reshape(-1, 1)), axis=1)\n",
        "\n",
        "df = pd.DataFrame(xy, columns=col_names + tgt_names)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "d4e2uj_OVnR_"
      },
      "outputs": [],
      "source": [
        "# 'Longitude', 'Latitude', 'HouseAge' - MinMaxScaler\n",
        "# 'Population', 'AveOccup', 'AveRooms', 'AveBedrms', 'MedInc' - PowerTransformer\n",
        "\n",
        "# standardize the target\n",
        "pipeline_min_max = Pipeline([\n",
        "    ('min_max_scaler', MinMaxScaler())\n",
        "])\n",
        "\n",
        "pipeline_power = Pipeline([\n",
        "    ('power_transformer', PowerTransformer())\n",
        "])\n",
        "\n",
        "\n",
        "pipeline_preproc = ColumnTransformer([\n",
        "    ('min_max_scaler', pipeline_min_max, ['Longitude', ]),\n",
        "    ('power_transformer', pipeline_power, ['Population', ]),\n",
        "],\n",
        "    remainder='passthrough')\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "CE4j_jJuVnR_"
      },
      "outputs": [],
      "source": [
        "pipeline_preproc.fit(df)\n",
        "df_preproc =  pd.DataFrame(pipeline_preproc.transform(df), columns=pipeline_preproc.get_feature_names_out())"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "QkTd5NZjVnR_"
      },
      "outputs": [],
      "source": [
        "df_preproc"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5CSi27XHVnR_"
      },
      "source": [
        "we can expand the pipeline to impute the missing values with the MICE method."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "JDiP3P9zVnR_"
      },
      "outputs": [],
      "source": [
        "pipeline_impute = Pipeline([\n",
        "    ('mice_imputer', mice)\n",
        "])\n",
        "\n",
        "pipeline_preproc = ColumnTransformer([\n",
        "    ('min_max_scaler', pipeline_min_max, ['Longitude', 'Latitude', 'HouseAge']),\n",
        "    ('power_transformer', pipeline_power, ['Population', 'AveOccup', 'AveRooms', 'AveBedrms', 'MedInc']),\n",
        "],\n",
        "    remainder='passthrough')\n",
        "\n",
        "pipeline = Pipeline([\n",
        "    ('preproc', pipeline_preproc),\n",
        "    ('impute', pipeline_impute),\n",
        "    ('pass', 'passthrough')\n",
        "])\n",
        "\n",
        "pipeline.fit(df)\n",
        "df_preproc =  pd.DataFrame(pipeline.transform(df), columns=pipeline.get_feature_names_out())\n",
        "df_preproc.head()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "hTsXYHQfVnR_"
      },
      "outputs": [],
      "source": [
        "# we need a custom function to make the pipeline work with the mice imputer function\n",
        "\n",
        "def mice_imputer(data_x):\n",
        "    if np.isnan(data_x).any():\n",
        "        return mice(data_x)\n",
        "    else:\n",
        "        return data_x\n",
        "\n",
        "pipeline_impute = Pipeline([\n",
        "    ('mice_imputer', FunctionTransformer(mice_imputer))\n",
        "])\n",
        "\n",
        "pipeline = Pipeline([\n",
        "    ('preproc', pipeline_preproc),\n",
        "    ('impute', pipeline_impute),\n",
        "    ('pass', 'passthrough')\n",
        "])\n",
        "\n",
        "pipeline.fit(df)\n",
        "df_preproc =  pd.DataFrame(pipeline.transform(df), columns=pipeline_preproc.get_feature_names_out())\n",
        "df_preproc.head()\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "A39t8SGOVnSA"
      },
      "outputs": [],
      "source": [
        "%%time\n",
        "# finally we can add a model to the pipeline. let's use a random forest regressor\n",
        "\n",
        "pipeline = Pipeline([\n",
        "    ('preproc', pipeline_preproc),\n",
        "    ('impute', pipeline_impute),\n",
        "    ('model', RandomForestRegressor())\n",
        "])\n",
        "\n",
        "tgt_name = tgt_names[0]\n",
        "df_train, df_test = train_test_split(df, test_size=0.2)\n",
        "pipeline.fit(df_train[col_names], df_train[tgt_name])\n",
        "\n",
        "y_pred = pipeline.predict(df_test[col_names])\n",
        "\n",
        "mse = mean_squared_error(df_test[tgt_name], y_pred)\n",
        "r2 = r2_score(df_test[tgt_name], y_pred)\n",
        "print(f'mse={mse:.3f}, r2={r2:.3f}')\n",
        "\n",
        "# plot the predicted vs true values\n",
        "\n",
        "fig, ax = plt.subplots(1, 1, figsize=(5, 5))\n",
        "sns.scatterplot(x=df_test[tgt_name], y=y_pred, ax=ax, alpha=0.3)\n",
        "ax.set_xlabel('True')\n",
        "ax.set_ylabel('Predicted')\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "NiGZdEOiVnSA"
      },
      "outputs": [],
      "source": [
        "%%time\n",
        "# we can also use the pipeline to do CV and grid search\n",
        "\n",
        "param_grid = {\n",
        "    'model__n_estimators': [10, 100, 200, 400],  # ToDo: test if 400 is ok on colab\n",
        "    'model__max_depth': [1, 20, 100],\n",
        "}\n",
        "\n",
        "grid_search = GridSearchCV(pipeline, param_grid, cv=2, n_jobs=2, verbose=2)\n",
        "\n",
        "grid_search.fit(df_train[col_names], df_train[tgt_name])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "-YuAeBTXVnSA"
      },
      "outputs": [],
      "source": [
        "print(grid_search.best_params_)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "tTg9tVYOVnSA"
      },
      "outputs": [],
      "source": [
        "grid_search.cv_results_"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "nuxizPHaVnSB"
      },
      "outputs": [],
      "source": [
        "res = grid_search.cv_results_\n",
        "cfgs = [str(p) for p in res['params']]\n",
        "scores = res['mean_test_score']\n",
        "scores_std = res['std_test_score']\n",
        "\n",
        "time = res['mean_fit_time']\n",
        "time_std = res['std_fit_time']\n",
        "\n",
        "# Results\n",
        "fig, ax = plt.subplots(1, 1, figsize=(10, 5))\n",
        "ax.errorbar(cfgs, scores, yerr=scores_std, fmt='o')\n",
        "ax.set_xticklabels(cfgs, rotation=90)\n",
        "ax.set_ylabel('mean test score (r2)')\n",
        "ax.set_xlabel('configurations')\n",
        "plt.show()\n",
        "\n",
        "\n",
        "# Time\n",
        "\n",
        "fig, ax = plt.subplots(1, 1, figsize=(10, 5))\n",
        "ax.errorbar(cfgs, time, yerr=time_std, fmt='o')\n",
        "ax.set_xticklabels(cfgs, rotation=90)\n",
        "ax.set_ylabel('mean fit time, sec')\n",
        "ax.set_xlabel('configurations')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "qWdWcyyUVnSB"
      },
      "outputs": [],
      "source": [
        "y_pred = grid_search.predict(df_test[col_names])\n",
        "\n",
        "mse = mean_squared_error(df_test[tgt_name], y_pred)\n",
        "r2 = r2_score(df_test[tgt_name], y_pred)\n",
        "\n",
        "print(f'mse={mse:.3f}, r2={r2:.3f}')\n",
        "\n",
        "# plot the predicted vs true values\n",
        "\n",
        "fig, ax = plt.subplots(1, 1, figsize=(5, 5))\n",
        "sns.scatterplot(x=df_test[tgt_name], y=y_pred, ax=ax, alpha=0.3)\n",
        "ax.set_xlabel('True')\n",
        "ax.set_ylabel('Predicted')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "GI2lYU2GVnSB"
      },
      "source": [
        "It is not advised to use the pipeline for the metric evaluation, e.g. see the discussion [here](https://stackoverflow.com/questions/43787107/use-a-metric-after-a-classifier-in-a-pipeline)    "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5AYAUOfsVnSB"
      },
      "source": [
        "# Step 4. Saving and Loading Models"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qjeYnC3vVnSB"
      },
      "source": [
        "The easiest is to use the `joblib` library, which is more efficient on objects that carry large numpy arrays internally as is often the case for fitted scikit-learn estimators.\n",
        "Saving the model with simple `pickle` can be problematic, as it can lead to compatibility issues between different Python versions."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "telTHZ8dVnSB"
      },
      "outputs": [],
      "source": [
        "# save the pipeline\n",
        "dump(pipeline, 'pipeline.joblib')\n",
        "\n",
        "# load the pipeline\n",
        "pipeline_loaded = load('pipeline.joblib')\n",
        "\n",
        "# use the loaded pipeline\n",
        "y_pred = pipeline_loaded.predict(df_test[col_names])\n",
        "\n",
        "mse = mean_squared_error(df_test[tgt_name], y_pred)\n",
        "print(f'mse={mse:.3f}')\n",
        "\n",
        "# plot the predicted vs true values\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "R_atz0VTVnSB"
      },
      "source": [
        "# Step 5. Data IO"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "c5OGQbGJVnSB"
      },
      "outputs": [],
      "source": [
        "# Save the data\n",
        "\n",
        "# csv is human-readable, but not the most efficient\n",
        "df.to_csv('california_housing.csv', index=False)\n",
        "\n",
        "# feather is the most efficient for pandas dataframes\n",
        "df.to_feather('california_housing.feather')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Ny4vd6IwVnSB"
      },
      "outputs": [],
      "source": [
        "# Load the data, testing the time\n",
        "\n",
        "%timeit df_csv = pd.read_csv('california_housing.csv')\n",
        "%timeit df_feather = pd.read_feather('california_housing.feather')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "1K-ka0wIVnSC"
      },
      "outputs": [],
      "source": [
        "# File sizes\n",
        "\n",
        "print(f'csv: {os.path.getsize(\"california_housing.csv\") / 1024 / 1024:.3f} MB')\n",
        "print(f'feather: {os.path.getsize(\"california_housing.feather\") / 1024 / 1024:.3f} MB')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tYrWZKfNVnSC"
      },
      "source": [
        "# Step 6. Data Visualization"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Gq-PL5r1VnSC"
      },
      "source": [
        "We will end this course with same message as we started it: the importance of data visualization.\n",
        "After all, the main goal of the data analysis is to make the data understandable and interpretable.\n",
        "Once we did all the ML work, we need to present the results in a way that is understandable - to the stakeholders, to the team, to the scientific community, to the public, and most importantly to ourselves.\n",
        "\n",
        "\n",
        "We can use plotly for interactive plots:\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "QOw4e9sQVnSC"
      },
      "outputs": [],
      "source": [
        "fig = px.scatter(df, x='AveRooms', y='AveBedrms', color='MedHouseVal', hover_data=['Latitude', 'Longitude'])"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "fig"
      ],
      "metadata": {
        "id": "bF2jeBV7j5eB"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "vkc3R1EZVnSC"
      },
      "outputs": [],
      "source": [
        "df['idx'] = df.index\n",
        "df['pred'] = grid_search.predict(df[col_names])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "v2rxz3f3VnSC",
        "scrolled": true
      },
      "outputs": [],
      "source": [
        "external_stylesheets = ['https://codepen.io/chriddyp/pen/bWLwgP.css']\n",
        "app = Dash('Houses', external_stylesheets=external_stylesheets)\n",
        "\n",
        "\n",
        "\n",
        "app.layout = html.Div([\n",
        "    html.Div([\n",
        "\n",
        "        html.Div([\n",
        "            dcc.Dropdown(\n",
        "                df.columns,\n",
        "                'AveRooms',\n",
        "                id='crossfilter-xaxis-column',\n",
        "            ),\n",
        "            dcc.RadioItems(\n",
        "                ['Linear', 'Log'],\n",
        "                'Linear',\n",
        "                id='crossfilter-xaxis-type',\n",
        "                labelStyle={'display': 'inline-block', 'marginTop': '5px'}\n",
        "            )\n",
        "        ],\n",
        "        style={'width': '32%', 'float': 'left', 'display': 'inline-block'}),\n",
        "\n",
        "        html.Div([\n",
        "            dcc.Dropdown(\n",
        "                df.columns,\n",
        "                'AveBedrms',\n",
        "                id='crossfilter-yaxis-column'\n",
        "            ),\n",
        "            dcc.RadioItems(\n",
        "                ['Linear', 'Log'],\n",
        "                'Linear',\n",
        "                id='crossfilter-yaxis-type',\n",
        "                labelStyle={'display': 'inline-block', 'marginTop': '5px'}\n",
        "            )\n",
        "        ], style={'width': '32%', 'float': 'left', 'display': 'inline-block'}),\n",
        "\n",
        "        html.Div([\n",
        "            dcc.Dropdown(\n",
        "                df.columns,\n",
        "                'MedHouseVal',\n",
        "                id='crossfilter-zaxis-column'\n",
        "            )\n",
        "        ], style={'width': '32%', 'float': 'left', 'display': 'inline-block'})\n",
        "    ], style={\n",
        "        'padding': '10px 5px'\n",
        "    }),\n",
        "\n",
        "    html.Div([\n",
        "        dcc.Graph(\n",
        "            id='crossfilter-indicator-scatter',\n",
        "            #hoverData={'points': [{'customdata': None}]},\n",
        "            selectedData={'points': [{'customdata': None}]}\n",
        "\n",
        "        )\n",
        "    ], style={'width': '49%', 'display': 'inline-block', 'padding': '0 20'}),\n",
        "    html.Div([\n",
        "        dcc.Graph(id='x-plot'),\n",
        "        dcc.Graph(id='y-plot'),\n",
        "    ], style={'display': 'inline-block', 'width': '49%'}),\n",
        "\n",
        "    html.Div(dcc.Slider(\n",
        "        df['MedHouseVal'].min(),\n",
        "        df['MedHouseVal'].max(),\n",
        "        step=None,\n",
        "        id='crossfilter-houseval--slider',\n",
        "        value=df['MedHouseVal'].max(),\n",
        "        marks={str(v): str(v) for v in range(int(df['MedHouseVal'].max()+1))}\n",
        "    ), style={'width': '49%', 'padding': '0px 20px 20px 20px'})\n",
        "])\n",
        "\n",
        "\n",
        "@callback(\n",
        "    Output('crossfilter-indicator-scatter', 'figure'),\n",
        "    Input('crossfilter-xaxis-column', 'value'),\n",
        "    Input('crossfilter-yaxis-column', 'value'),\n",
        "    Input('crossfilter-zaxis-column', 'value'),\n",
        "    Input('crossfilter-xaxis-type', 'value'),\n",
        "    Input('crossfilter-yaxis-type', 'value'),\n",
        "    Input('crossfilter-houseval--slider', 'value'))\n",
        "def update_graph(xaxis_column_name, yaxis_column_name, zaxis_column_name,\n",
        "                 xaxis_type, yaxis_type,\n",
        "                 price_value):\n",
        "    dff = df[df['MedHouseVal'] < price_value]\n",
        "\n",
        "\n",
        "\n",
        "    fig = px.scatter(x=dff[xaxis_column_name],\n",
        "                     y=dff[yaxis_column_name],\n",
        "                     color=dff[zaxis_column_name],\n",
        "                     hover_name=dff['MedHouseVal']\n",
        "                     )\n",
        "\n",
        "    fig.update_traces(customdata=dff['idx'])\n",
        "\n",
        "    fig.update_xaxes(title=xaxis_column_name, type='linear' if xaxis_type == 'Linear' else 'log')\n",
        "\n",
        "    fig.update_yaxes(title=yaxis_column_name, type='linear' if yaxis_type == 'Linear' else 'log')\n",
        "\n",
        "    fig.update_layout(margin={'l': 40, 'b': 40, 't': 10, 'r': 0}, hovermode='closest')\n",
        "\n",
        "    return fig\n",
        "\n",
        "\n",
        "def create_plot(dff, axis_type, axis):\n",
        "\n",
        "    fig = px.scatter(dff, x='MedHouseVal', y=axis)\n",
        "\n",
        "    fig.update_traces(mode='markers')\n",
        "\n",
        "    fig.update_xaxes(showgrid=False)\n",
        "\n",
        "    fig.update_yaxes(type='linear' if axis_type == 'Linear' else 'log')\n",
        "\n",
        "    fig.add_annotation(x=0, y=0.85, xanchor='left', yanchor='bottom',\n",
        "                       xref='paper', yref='paper', showarrow=False, align='left',\n",
        "                       text=axis)\n",
        "\n",
        "    fig.update_layout(height=225, margin={'l': 20, 'b': 30, 'r': 10, 't': 10})\n",
        "\n",
        "    return fig\n",
        "\n",
        "\n",
        "@callback(\n",
        "    Output('x-plot', 'figure'),\n",
        "    Input('crossfilter-indicator-scatter', 'selectedData'),\n",
        "    Input('crossfilter-xaxis-column', 'value'),\n",
        "    Input('crossfilter-xaxis-type', 'value'))\n",
        "def update_x_timeseries(hoverData, xaxis_column_name, axis_type):\n",
        "    idx = [p['customdata'] for p in hoverData['points']]\n",
        "    idx = [i for i in idx if i is not None]\n",
        "    dff = df.iloc[idx]\n",
        "    return create_plot(dff, axis_type, xaxis_column_name)\n",
        "\n",
        "\n",
        "@callback(\n",
        "    Output('y-plot', 'figure'),\n",
        "    Input('crossfilter-indicator-scatter', 'selectedData'),\n",
        "    Input('crossfilter-yaxis-column', 'value'),\n",
        "    Input('crossfilter-yaxis-type', 'value'))\n",
        "def update_y_timeseries(hoverData, yaxis_column_name, axis_type):\n",
        "    idx = [p['customdata'] for p in hoverData['points']]\n",
        "    idx = [i for i in idx if i is not None]\n",
        "    dff = df.iloc[idx]\n",
        "    return create_plot(dff, axis_type, yaxis_column_name)\n",
        "\n",
        "\n",
        "app.run(debug=True)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "lF17p0bQVnSD"
      },
      "source": [
        "# Workshop"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "KzDpaIGRVnSD"
      },
      "source": [
        "Use your data to exercise the same steps as in the notebook.\n",
        "\n",
        "Alternatively, you can use the extended house price data: Subset of the Ames Houses dataset: http://jse.amstat.org/v19n3/decock.pdf\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "45F-ShhqVnSD"
      },
      "outputs": [],
      "source": [
        "# load the data\n",
        "path = os.path.abspath('.')+'/colab_material.tgz'\n",
        "url = 'https://github.com/neworldemancer/DSF5/raw/master/colab_material.tgz'\n",
        "\n",
        "r = requests.get(url)\n",
        "\n",
        "with open(path, 'wb') as f:\n",
        "    f.write(r.content)\n",
        "\n",
        "tar = tarfile.open(path, \"r:gz\")\n",
        "tar.extractall()\n",
        "tar.close()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "H2jaGbJyVnSD"
      },
      "outputs": [],
      "source": [
        "path = 'data/AmesHousing.csv'\n",
        "df = pd.read_csv(path, na_values=('NaN', ''), keep_default_na=False)\n",
        "rename_dict = {k:k.replace(' ', '').replace('/', '') for k in df.keys()}  # simplify the column names\n",
        "df.rename(columns=rename_dict, inplace=True)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "88FoCCXxVnSD"
      },
      "outputs": [],
      "source": [
        "df.head()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "14xUgUpyVnSD"
      },
      "source": [
        "# Further reading"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "9R4ZaD3rVnSD"
      },
      "source": [
        "* [scikit-learn](https://scikit-learn.org/stable/)\n",
        "* [Clustering with visualizations](https://dashee87.github.io/data%20science/general/Clustering-with-Scikit-with-GIFs)\n",
        "* [hierarchical clustering](https://scikit-learn.org/stable/modules/clustering.html#hierarchical-clustering)\n",
        "* [dash](https://dash.plotly.com/)\n",
        "* [plotly](https://plotly.com/python/)\n",
        "* [streamlit](https://streamlit.io/)\n",
        "* [MLFlow](https://mlflow.org/)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "RY1ZTco4VnSE"
      },
      "source": [
        "# Courses and support at DSL"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": false,
        "id": "NvbCl8R0VnSE"
      },
      "source": [
        "* [Data Science Lab](https://www.dsl.unibe.ch)\n",
        "* [Data Science Lab Walk-in](https://www.dsl.unibe.ch/support/walk_in/)"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "provenance": [],
      "toc_visible": true
    },
    "kernelspec": {
      "display_name": "Python 3 (ipykernel)",
      "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.10.12"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}