docs/source/examples/data_drift/data_drift_examples.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data Drift Detector Examples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook contains examples that show how to set up, run, and produce output from detectors in the \n",
    "data_drift module. The parameters aren't necessarily tuned for best \n",
    "performance for the input data, just notional. These detectors are generally\n",
    "applied to the whole feature set for a given data source. \n",
    "\n",
    "The example data for PCA-CD and KDQ-Tree, Circle, is a synthetic data source, where drift\n",
    "occurs in both var1, var2, and the conditional distributions P(y|var1) and\n",
    "P(y|var2). The drift occurs from index 1000 to 1250, and affects 66% of the\n",
    "sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Imports ##\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import plotly.express as px\n",
    "from menelaus.data_drift.cdbd import CDBD\n",
    "from menelaus.data_drift.hdddm import HDDDM\n",
    "from menelaus.data_drift import PCACD\n",
    "from menelaus.data_drift import KdqTreeStreaming, KdqTreeBatch\n",
    "from menelaus.data_drift import NNDVI\n",
    "from menelaus.datasets import make_example_batch_data, fetch_circle_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Import Data ##\n",
    "\n",
    "# Import CDBD and HDDDM Data\n",
    "example_data = make_example_batch_data()\n",
    "\n",
    "# Import PCA-CD and KDQ-Tree Data (circle dataset)\n",
    "circle_data = fetch_circle_data()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Confidence Distribution Batch Detection (CDBD)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This section details how to setup, run, and produce plots for CDBD. This script\n",
    "monitors the feature \"confidence\", simulated confidence scores output by a\n",
    "classifier. Drift occurs in 2018 and persists through 2021. See\n",
    "`make_example_batch_data` for more info.\n",
    "\n",
    "CDBD must be setup and run with batches of data containing 1 variable. \n",
    "\n",
    "Plots include: \n",
    "\n",
    "- A line plot visualizing test statistics for detection of drift"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Setup ##\n",
    "\n",
    "# Set up reference and test batches, using 2007 as reference year\n",
    "reference = pd.DataFrame(example_data[example_data.year == 2007].loc[:, \"confidence\"])\n",
    "all_test = example_data[example_data.year != 2007]\n",
    "\n",
    "# Run CDBD\n",
    "cdbd = CDBD(subsets=8)\n",
    "cdbd.set_reference(reference)\n",
    "\n",
    "# Store drift for test statistic plot\n",
    "detected_drift = []\n",
    "for year, subset_data in all_test.groupby(\"year\"):\n",
    "    cdbd.update(pd.DataFrame(subset_data.loc[:, \"confidence\"]))\n",
    "    detected_drift.append(cdbd.drift_state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Drift occurs in 2018 and persists through end of dataset. CDBD identifies drift\n",
    "occurring in 2019, one year late. It alerts to a false alarm in 2012.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Plot Line Graph ##\n",
    "# Calculate divergences for all years in dataset\n",
    "years = list(example_data.year.value_counts().index[1:])\n",
    "kl_divergence = [\n",
    "    ep - th for ep, th in zip(cdbd.epsilon_values.values(), cdbd.thresholds.values())\n",
    "]\n",
    "\n",
    "# Remove potential infs that arise because of small confidence scores\n",
    "kl_divergence = [\n",
    "    x if np.isnan(x) == False and np.isinf(x) == False else 0 for x in kl_divergence\n",
    "]\n",
    "\n",
    "# Plot KL Divergence against Year, along with detected drift\n",
    "plot_data = pd.DataFrame(\n",
    "    {\"Year\": years, \"KL Divergence\": kl_divergence, \"Detected Drift\": detected_drift}\n",
    ")\n",
    "sns.set_style(\"white\")\n",
    "plt.figure(figsize=(20, 6))\n",
    "plt.plot(\"Year\", \"KL Divergence\", data=plot_data, label=\"KL Divergence\", marker=\".\")\n",
    "plt.grid(False, axis=\"x\")\n",
    "plt.xticks(years, fontsize=16)\n",
    "plt.yticks(fontsize=16)\n",
    "plt.title(\"CDBD Test Statistics\", fontsize=22)\n",
    "plt.ylabel(\"KL Divergence\", fontsize=18)\n",
    "plt.xlabel(\"Year\", fontsize=18)\n",
    "plt.ylim([min(kl_divergence) - 0.02, max(kl_divergence) + 0.02])\n",
    "for _, t in enumerate(plot_data.loc[plot_data[\"Detected Drift\"] == \"drift\"][\"Year\"]):\n",
    "    plt.axvspan(\n",
    "        t - 0.2, t + 0.2, alpha=0.5, color=\"red\", label=(\"Drift Detected\" if _ == 0 else None)\n",
    "    )\n",
    "plt.legend()\n",
    "plt.axhline(y=0, color=\"orange\", linestyle=\"dashed\")\n",
    "\n",
    "plt.show()\n",
    "# plt.savefig(\"example_CDBD_test_statistics.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Custom Divergence Metric ### \n",
    "\n",
    "# Define divergence function\n",
    "def distance_metric(reference_histogram, test_histogram):\n",
    "\n",
    "    # Convert inputs to appropriate datatype \n",
    "    ref = np.array(reference_histogram[0])\n",
    "    test = np.array(test_histogram[0])\n",
    "\n",
    "    return np.sqrt(np.sum(np.square(ref-test)))\n",
    "\n",
    "# Test self-defined divergence metric \n",
    "cdbd = CDBD(\n",
    "    divergence=distance_metric,\n",
    "    detect_batch=1,\n",
    "    statistic=\"stdev\",\n",
    "    significance=0.05,\n",
    "    subsets=5,\n",
    ")\n",
    "\n",
    "cdbd.set_reference(reference)\n",
    "cdbd.update(pd.DataFrame(example_data[example_data.year == 2008].loc[:, \"confidence\"]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hellinger Distance Drift Detection Method (HDDDM)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This section details how to setup, run, and produce plots for HDDDM, using both\n",
    "numeric and categorical data. Drift occurs in 2009, 2012, 2015, 2018, and 2021.\n",
    "Drift in 2018 persists through 2021. See `make_example_batch_data` for more info.\n",
    "\n",
    "It additionally contains an example of using a custom divergence function. \n",
    "\n",
    "HDDDM must be setup and run with batches of data. \n",
    "\n",
    "Plots include:\n",
    "\n",
    "- A line plot visualizing test statistics for detection of drift \n",
    "\n",
    "- A heatmap visualizing \"where\" drift is occuring, showing features for each\n",
    "  year's test batch with the greatest Hellinger distance from the reference\n",
    "  batch."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Setup ##\n",
    "\n",
    "# Set up reference and test batches, using 2007 as reference year\n",
    "# -2 indexing removes columns \"drift\" and \"confidence\"\n",
    "reference = example_data[example_data.year == 2007].iloc[:, 1:-2]\n",
    "all_test = example_data[example_data.year != 2007]\n",
    "\n",
    "#  Setup HDDDM\n",
    "np.random.seed(1)\n",
    "hdddm = HDDDM(subsets=8)\n",
    "\n",
    "# Store epsilons per feature for heatmap\n",
    "years = all_test.year.unique()\n",
    "heatmap_data = pd.DataFrame(columns = years)\n",
    "\n",
    "# Store drift for test statistic plot\n",
    "detected_drift = []\n",
    "\n",
    "# Run HDDDM\n",
    "hdddm.set_reference(reference)\n",
    "for year, subset_data in example_data[example_data.year != 2007].groupby(\"year\"):\n",
    "    hdddm.update(subset_data.iloc[:, 1:-2])\n",
    "    heatmap_data[year] = hdddm.feature_epsilons\n",
    "    detected_drift.append(hdddm.drift_state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "HDDDM identifies drifts in 2009, 2010, 2012, 2019, 2021. These drifts involve\n",
    "a change in mean or variance. Drift in 2010 is likely identified as the\n",
    "distribution returns to state prior to 2009 drift. Drift in 2015, a change in\n",
    "correlation, is undetected. Drift in 2018 is detected one year late.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Plot Line Graph ##\n",
    "h_distances = [\n",
    "    ep - th for ep, th in zip(hdddm.epsilon_values.values(), hdddm.thresholds.values())\n",
    "]\n",
    "\n",
    "# Plot Hellinger Distance against Year, along with detected drift\n",
    "plot_data = pd.DataFrame(\n",
    "    {\"Year\": years, \"Hellinger Distance\": h_distances, \"Detected Drift\": detected_drift}\n",
    ")\n",
    "sns.set_style(\"white\")\n",
    "plt.figure(figsize=(20, 6))\n",
    "plt.plot(\n",
    "    \"Year\", \"Hellinger Distance\", data=plot_data, label=\"Hellinger Distance\", marker=\".\"\n",
    ")\n",
    "plt.grid(False, axis=\"x\")\n",
    "plt.xticks(years, fontsize=16)\n",
    "plt.yticks(fontsize=16)\n",
    "plt.title(\"HDDDM Test Statistics\", fontsize=22)\n",
    "plt.ylabel(\"Hellinger Distance\", fontsize=18)\n",
    "plt.xlabel(\"Year\", fontsize=18)\n",
    "plt.ylim([min(h_distances) - 0.02, max(h_distances) + 0.02])\n",
    "for _, t in enumerate(plot_data.loc[plot_data[\"Detected Drift\"] == \"drift\"][\"Year\"]):\n",
    "    plt.axvspan(\n",
    "        t - 0.2, t + 0.2, alpha=0.5, color=\"red\", label=(\"Drift Detected\" if _ == 0 else None)\n",
    "    )\n",
    "plt.legend()\n",
    "plt.axhline(y=0, color=\"orange\", linestyle=\"dashed\")\n",
    "\n",
    "plt.show()\n",
    "# plt.savefig(\"example_HDDDM_test_statistics.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Drift in feature B is detected in 2009 and 2010 (as it reverts to normal).\n",
    "\n",
    "- Drift in feature D is detected in 2012 and 2013 (as it reverts to normal).\n",
    "\n",
    "- Drift in feature H is detected in 2019. Drift in feature J is detected in 2021.\n",
    "\n",
    "- The undetected drift occurs in 2015 in the correlations between features E and F.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Plot Heatmap ##\n",
    "sns.set_style(\"whitegrid\")\n",
    "sns.set(rc={\"figure.figsize\": (15, 8)})\n",
    "\n",
    "# Setup plot\n",
    "\n",
    "# Setup plot\n",
    "grid_kws = {\"height_ratios\": (0.9, 0.05), \"hspace\": 0.3}\n",
    "f, (ax, cbar_ax) = plt.subplots(2, gridspec_kw=grid_kws)\n",
    "coloring = sns.cubehelix_palette(start=0.8, rot=-0.5, as_cmap=True)\n",
    "ax = sns.heatmap(\n",
    "    heatmap_data,\n",
    "    ax=ax,\n",
    "    cmap=coloring,\n",
    "    xticklabels=heatmap_data.columns,\n",
    "    yticklabels=heatmap_data.index,\n",
    "    linewidths=0.5,\n",
    "    cbar_ax=cbar_ax,\n",
    "    cbar_kws={\"orientation\": \"horizontal\"},\n",
    ")\n",
    "ax.set_title('HDDDM Feature Heatmap')\n",
    "ax.set(xlabel=\"Years\", ylabel=\"Features\")\n",
    "ax.collections[0].colorbar.set_label(\"Difference in Hellinger Distance\")\n",
    "ax.set_yticklabels(ax.get_yticklabels(), rotation=0)\n",
    "\n",
    "plt.show()\n",
    "# plt.savefig(\"example_HDDDM_feature_heatmap.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Custom Divergence Metric ### \n",
    "\n",
    "# Define divergence function\n",
    "def distance_metric(reference_histogram, test_histogram):\n",
    "\n",
    "    # Convert inputs to appropriate datatype \n",
    "    ref = np.array(reference_histogram[0])\n",
    "    test = np.array(test_histogram[0])\n",
    "\n",
    "    return np.sqrt(np.sum(np.square(ref-test)))\n",
    "\n",
    "# Test self-defined divergence metric \n",
    "hdddm = HDDDM(\n",
    "    divergence=distance_metric,\n",
    "    detect_batch=1,\n",
    "    statistic=\"stdev\",\n",
    "    significance=0.05,\n",
    "    subsets=5,\n",
    ")\n",
    "\n",
    "hdddm.set_reference(reference)\n",
    "hdddm.update(example_data[example_data.year == 2008].iloc[:, 1:-2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PCA-Based Change Detection (PCA-CD)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PCA-CD is a drift detector that transforms the passed data into its principal components, then watches the transformed data for signs of drift by monitoring the KL-divergence via the Page-Hinkley algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Setup ##\n",
    "\n",
    "pca_cd = PCACD(window_size=50, divergence_metric=\"intersection\")\n",
    "\n",
    "# set up dataframe to record results\n",
    "status = pd.DataFrame(columns=[\"index\", \"var1\", \"var2\", \"drift_detected\"])\n",
    "\n",
    "# Put together a dataframe of several features, each of which abruptly changes\n",
    "# at index 1000.\n",
    "np.random.seed(1)\n",
    "size = 1000\n",
    "data = pd.concat(\n",
    "    [\n",
    "        pd.DataFrame(\n",
    "            [\n",
    "                np.random.normal(1, 10, size),\n",
    "                np.random.uniform(1, 2, size),\n",
    "                np.random.normal(0, 1, size),\n",
    "            ]\n",
    "        ).T,\n",
    "        pd.DataFrame(\n",
    "            [\n",
    "                np.random.normal(9, 10, size),\n",
    "                np.random.normal(1, 3, size),\n",
    "                np.random.gamma(20, 30, size),\n",
    "            ]\n",
    "        ).T,\n",
    "    ]\n",
    ")\n",
    "\n",
    "# Update the drift detector with each new sample\n",
    "for i in range(len(circle_data)):\n",
    "    pca_cd.update(data.iloc[[i]])\n",
    "    status.loc[i] = [i, data.iloc[i, 0], data.iloc[i, 1], pca_cd.drift_state]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Plotting ##\n",
    "\n",
    "# Plot the features and the drift\n",
    "plt.figure(figsize=(20, 6))\n",
    "plt.scatter(status.index, status.var2, label=\"Var 2\")\n",
    "plt.scatter(status.index, status.var1, label=\"Var 1\", alpha=0.5)\n",
    "plt.grid(False, axis=\"x\")\n",
    "plt.xticks(fontsize=16)\n",
    "plt.yticks(fontsize=16)\n",
    "plt.title(\"PCA-CD Test Results\", fontsize=22)\n",
    "plt.ylabel(\"Value\", fontsize=18)\n",
    "plt.xlabel(\"Index\", fontsize=18)\n",
    "ylims = min(status.var1.min(), status.var2.min()), max(\n",
    "    status.var1.max(), status.var1.max()\n",
    ")\n",
    "plt.ylim(ylims)\n",
    "\n",
    "# Draw red lines that indicate where drift was detected\n",
    "plt.vlines(\n",
    "    x=status.loc[status[\"drift_detected\"] == \"drift\"][\"index\"],\n",
    "    ymin=ylims[0],\n",
    "    ymax=ylims[1],\n",
    "    label=\"Drift Detected\",\n",
    "    color=\"red\",\n",
    ")\n",
    "\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PCA-CD detects this very abrupt drift within a few samples of its induction.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.show()\n",
    "# plt.savefig(\"example_PCA_CD.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## KDQ-Tree Detection Method (Streaming Setting)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "KdqTree monitors incoming features by constructing a tree which partitions the feature-space, and then monitoring a divergence statistic that is defined over that partition. It watches data within a sliding window of a particular size. When that window is full, it builds the reference tree. As the window moves forward, point-by-point, the data in that new window is compared against the reference tree to detect drift."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Setup ##\n",
    "\n",
    "# kdqTree does use bootstrapping to define its critical thresholds, so setting\n",
    "# the seed is important to reproduce exact behavior.\n",
    "np.random.seed(1)\n",
    "\n",
    "# Note that the default input_type for KDQTree is \"stream\".\n",
    "# The window size, corresponding to the portion of the stream which KDQTree\n",
    "# monitors, must be specified.\n",
    "det = KdqTreeStreaming(window_size=500, alpha=0.05, bootstrap_samples=500, count_ubound=50)\n",
    "\n",
    "# setup DF to record results\n",
    "status = pd.DataFrame(columns=[\"index\", \"var1\", \"var2\", \"drift_detected\"])\n",
    "\n",
    "# iterate through X data and run detector\n",
    "data = circle_data[[\"var1\", \"var2\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Plotting ##\n",
    "\n",
    "plot_data = {}\n",
    "for i in range(len(circle_data)):\n",
    "    det.update(data.iloc[[i]])\n",
    "    status.loc[i] = [i, data.iloc[i, 0], data.iloc[i, 1], det.drift_state]\n",
    "    if det.drift_state is not None:\n",
    "        # capture the visualization data\n",
    "        plot_data[i] = det.to_plotly_dataframe()\n",
    "\n",
    "plt.figure(figsize=(20, 6))\n",
    "plt.scatter(\"index\", \"var2\", data=status, label=\"var2\")\n",
    "plt.scatter(\"index\", \"var1\", data=status, label=\"var1\", alpha=0.5)\n",
    "plt.grid(False, axis=\"x\")\n",
    "plt.xticks(fontsize=16)\n",
    "plt.yticks(fontsize=16)\n",
    "plt.title(\"KDQ Tree Test Results\", fontsize=22)\n",
    "plt.ylabel(\"Value\", fontsize=18)\n",
    "plt.xlabel(\"Index\", fontsize=18)\n",
    "ylims = [-0.05, 1.05]\n",
    "plt.ylim(ylims)\n",
    "\n",
    "drift_start, drift_end = 1000, 1250\n",
    "plt.axvspan(drift_start, drift_end, alpha=0.5, label=\"Drift Induction Window\")\n",
    "\n",
    "# Draw red lines that indicate where drift was detected\n",
    "plt.vlines(\n",
    "    x=status.loc[status[\"drift_detected\"] == \"drift\"][\"index\"],\n",
    "    ymin=ylims[0],\n",
    "    ymax=ylims[1],\n",
    "    label=\"Drift Detected\",\n",
    "    color=\"red\",\n",
    ")\n",
    "\n",
    "plt.legend()\n",
    "plt.show()\n",
    "# plt.savefig(\"example_streaming_kdqtree_feature_stream.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given a window_size of 500, with only the two input features, KdqTree detects\n",
    "a change after about half of the data within its window is in the new regime.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we save off the to_plotly_dataframe at each drift detection, we can display\n",
    "the Kulldorff Spatial Scan Statistic (KSS) for each. Higher values of KSS\n",
    "indicate that a given region of the data space has greater divergence between\n",
    "the reference and test data.\n",
    "\n",
    "Note that the structure of the particular tree depends on the reference data\n",
    "and the order of the columns within the dataframe!\n",
    "\n",
    "Since this data only contains two features, the tree is relatively\n",
    "shallow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Kulldorff Spatial Scan Statistic (KSS) ##\n",
    "for title, df_plot in plot_data.items():\n",
    "    fig = px.treemap(\n",
    "        data_frame=df_plot,\n",
    "        names=\"name\",\n",
    "        ids=\"idx\",\n",
    "        parents=\"parent_idx\",\n",
    "        color=\"kss\",\n",
    "        color_continuous_scale=\"blues\",\n",
    "        title=f\"Index {title}\",\n",
    "    )\n",
    "    fig.update_traces(root_color=\"lightgrey\")\n",
    "    fig.show()\n",
    "    # fig.write_html(f\"example_streaming_kdqtree_treemap_{title}.html\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## KDQ-Tree Detection Method (Batch Setting)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example shows up how to set up, run, and produce output from the kdq-Tree\n",
    "detector, specifically in the batch data setting. The parameters aren't \n",
    "necessarily tuned for best performance, just notional.\n",
    "\n",
    "Drift in the example dataset occurs in 2009, 2012, 2015, 2018, and 2021. Drift\n",
    "in 2018 persists through 2021. See\n",
    "`make_example_batch_data` for more details.\n",
    "\n",
    "This example takes roughly a minute to run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Setup ##\n",
    "\n",
    "# kdq-Tree does use bootstrapping to define its critical thresholds, so setting\n",
    "# the seed is important to reproduce exact behavior.\n",
    "np.random.seed(123)\n",
    "\n",
    "df = example_data\n",
    "\n",
    "# Capture the column which tells us when drift truly occurred\n",
    "drift_years = df.groupby(\"year\")[\"drift\"].apply(lambda x: x.unique()[0]).reset_index()\n",
    "\n",
    "# Because the drift in 2009, 2012, and 2016 is intermittent - it reverts\n",
    "# back to the prior distribution in the subsequent year - we should also detect\n",
    "# drift in 2010, 2013, and 2016. So:\n",
    "drift_years.loc[drift_years[\"year\"].isin([2010, 2013, 2016]), \"drift\"] = True\n",
    "\n",
    "\n",
    "df.drop(columns=[\"cat\", \"confidence\", \"drift\"], inplace=True)\n",
    "\n",
    "\n",
    "plot_data = {}\n",
    "status = pd.DataFrame(columns=[\"year\", \"drift\"])\n",
    "det = KdqTreeBatch()\n",
    "\n",
    "# Set up reference batch, using 2007 as reference year\n",
    "det.set_reference(df[df.year == 2007].drop(columns=['year']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Batch the data by year and run kdqTree\n",
    "for group, sub_df in df[df.year != 2007].groupby(\"year\"):\n",
    "    det.update(sub_df.drop(columns=[\"year\"]))\n",
    "    status = pd.concat(\n",
    "        [status, pd.DataFrame({\"year\": [group], \"drift\": [det.drift_state]})],\n",
    "        axis=0,\n",
    "        ignore_index=True,\n",
    "    )\n",
    "    if det.drift_state is not None:\n",
    "        # capture the visualization data\n",
    "        plot_data[group] = det.to_plotly_dataframe()\n",
    "        \n",
    "        # option to specify reference batch to be any year \n",
    "        #det.set_reference(df[df.year == XXXX])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print out the true drift status, and that according to the detector.\n",
    "The detector successfully identifies drift in every year but 2018;\n",
    "that's the first year where persistent drift, from 2018-2021, was induced.\n",
    "The detector picks it up in 2019, the second year of persistent drift."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(\n",
    "    status.merge(drift_years, how=\"left\", on=\"year\", suffixes=[\"_kdqTree\", \"_true\"])\n",
    "    .replace({True: \"drift\", False: None})\n",
    "    # .to_csv(\"example_kdqtree_drift_comparison.csv\", index=False)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we save off the to_plotly_dataframe at each drift detection, we can display\n",
    "the Kulldorff Spatial Scan Statistic (KSS) for each. Higher values of KSS\n",
    "indicate that a given region of the data space has greater divergence between\n",
    "the reference and test data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for year, df_plot in plot_data.items():\n",
    "    fig = px.treemap(\n",
    "        data_frame=df_plot,\n",
    "        names=\"name\",\n",
    "        ids=\"idx\",\n",
    "        parents=\"parent_idx\",\n",
    "        color=\"kss\",\n",
    "        color_continuous_scale=\"blues\",\n",
    "        title=year,\n",
    "    )\n",
    "    fig.update_traces(root_color=\"lightgrey\")\n",
    "    fig.show()\n",
    "    # fig.write_html(f\"example_kdqtree_treemap_{year}.html\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "We can see that the regions of greatest drift do line up with at least one of\n",
    "the items that were modified in a given year.\n",
    "\n",
    "For reference, the detailed descriptions of drift induction:\n",
    "\n",
    "- Drift 1: change the mean & var of item B in 2009, means will revert for 2010 on\n",
    "\n",
    "- Drift 2: change the variance of item c and d in 2012 by replacing some with the mean\n",
    "    keep same mean as other years, revert by 2013\n",
    "\n",
    "- Drift 3: change the correlation of item e and f in 2015 (go from correlation of 0 to correlation of 0.5)\n",
    "\n",
    "- Drift 4: change mean and var of H and persist it from 2018 on\n",
    "\n",
    "- Drift 5: change mean and var just for a year of J in 2021\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Nearest-Neighbor Density Variation Identification (NN-DVI)\n",
    "\n",
    "This example shows up how to set up, run, and produce output from the NN-DVI\n",
    "detector, which is implemented to run in the batch setting. NN-DVI uses k-nearest \n",
    "neighbors to partition the data, then uses measures defined with the adjacency \n",
    "matrix for the reference and test data to determine whether drift has occurred."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = make_example_batch_data()\n",
    "drift_years = data.groupby(\"year\")[\"drift\"].apply(lambda x: x.unique()[0]).reset_index()\n",
    "drift_years.loc[drift_years[\"year\"].isin([2010, 2013, 2016]), \"drift\"] = True\n",
    "\n",
    "data_grouped = data.groupby('year')\n",
    "\n",
    "# Note that a small subsample is used, to decrease the runtime for the purpose\n",
    "# of example.\n",
    "batches = {year: (group.sample(frac=.1)\n",
    "\t\t\t\t\t\t.drop(['year', 'cat', 'confidence', 'drift'], axis=1)\n",
    "\t\t\t\t\t\t.values)\n",
    "\t\t   for year, group in data_grouped}\n",
    "\n",
    "status = pd.DataFrame(columns=[\"year\", \"drift\"])\n",
    "\n",
    "det = NNDVI(k_nn=2, sampling_times=50)\n",
    "det.set_reference(batches.pop(2007))\n",
    "\n",
    "for year, batch in batches.items():\n",
    "\tdet.update(batch)\n",
    "\tstatus = pd.concat([status, pd.DataFrame({\"year\":[year], \"drift\":[det.drift_state]})], \n",
    "\t\t\t\t\t\taxis=0, \n",
    "                        ignore_index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using this small subsample of data runs quickly, but yields inaccurate results \n",
    "for the test data. A run with the full data takes much longer, but gives the \n",
    "correct result in 9 of 13 years."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(\n",
    "    status.merge(drift_years, how=\"left\", on=\"year\", suffixes=[\"_nndvi\", \"_true\"])\n",
    "    .replace({True: \"drift\", False: None})\n",
    "    # .to_csv(\"example_nndvi_drift_comparison.csv\", index=False)\n",
    ")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.2 ('molten_env')",
   "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.2"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "cebda30e5482c3c83cb9c7c7d8ef1c1d1f67dcf93284b9e88f69b68f560a7bf6"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}