{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "r2d-title",
   "metadata": {},
   "source": [
    "# Rare-event identification on a 2D field\n",
    "\n",
    "A minimal, self-contained rare-event identification (REI) run: build a small 2D\n",
    "von-Mises stress field with a few high-stress patches, cluster it with\n",
    "`ClusterAnalysisIndicator`, and visualize the field next to the recovered\n",
    "clusters. This tutorial is pure Python and runs from `pip install graintrace`.\n",
    "\n",
    "See the algorithm page :doc:`/algorithms/rare-event-identification` and the API\n",
    ":class:`~graintrace.ClusterAnalysisIndicator`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "r2d-imports",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-09T23:29:59.631975Z",
     "iopub.status.busy": "2026-10-09T23:29:59.631501Z",
     "iopub.status.idle": "2026-10-09T23:30:00.232919Z",
     "shell.execute_reply": "2026-10-09T23:30:00.231328Z"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colors import ListedColormap, BoundaryNorm\n",
    "\n",
    "from graintrace.cluster_indicator import ClusterAnalysisIndicator\n",
    "from graintrace.similarity_metric_library import SimilarityMetricLibrary"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "r2d-gen-md",
   "metadata": {},
   "source": [
    "## Build a synthetic field\n",
    "\n",
    "Each grid point gets a Cauchy-stress row (`sxx..syz`). The background is low\n",
    "stress; a handful of circular patches are raised to a high value, so the rare\n",
    "regions are known ahead of time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "r2d-gen",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-09T23:30:00.237986Z",
     "iopub.status.busy": "2026-10-09T23:30:00.237717Z",
     "iopub.status.idle": "2026-10-09T23:30:00.265298Z",
     "shell.execute_reply": "2026-10-09T23:30:00.263631Z"
    }
   },
   "outputs": [],
   "source": [
    "def make_vms_grid(path, nx=30, ny=30, n_patches=10, vm_low=50.0, vm_high=200.0,\n",
    "                  radius_range=(2, 5), seed=42):\n",
    "    rng = np.random.default_rng(seed)\n",
    "    vm = np.full((ny, nx), vm_low)\n",
    "    for _ in range(n_patches):\n",
    "        cx, cy = rng.integers(0, nx), rng.integers(0, ny)\n",
    "        r = rng.integers(*radius_range)\n",
    "        yy, xx = np.ogrid[:ny, :nx]\n",
    "        vm[(xx - cx) ** 2 + (yy - cy) ** 2 <= r * r] = vm_high\n",
    "    rows = []\n",
    "    eid = 0\n",
    "    for j in range(ny):\n",
    "        for i in range(nx):\n",
    "            eid += 1\n",
    "            t = vm[j, i]\n",
    "            rows.append(dict(id=eid, x=float(i), y=float(j), z=0.0,\n",
    "                             sxx=t + rng.normal(0, t * 0.05),\n",
    "                             syy=rng.normal(0, t * 0.02), szz=rng.normal(0, t * 0.02),\n",
    "                             sxy=rng.normal(0, t * 0.02), sxz=rng.normal(0, t * 0.02),\n",
    "                             syz=rng.normal(0, t * 0.02)))\n",
    "    pd.DataFrame(rows).to_csv(path, index=False)\n",
    "    return nx, ny\n",
    "\n",
    "csv_path = \"rei2d_field.csv\"\n",
    "nx, ny = make_vms_grid(csv_path)\n",
    "print(f\"wrote {csv_path}: {nx}x{ny} grid\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "r2d-cluster-md",
   "metadata": {},
   "source": [
    "## Cluster with a von-Mises similarity metric\n",
    "\n",
    "`ClusterAnalysisIndicator` runs a single-stage hierarchical clustering over a\n",
    "feature from a `SimilarityMetricLibrary` spec (here the von-Mises stress).\n",
    "`threshold` sets the distance cut of the dendrogram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "r2d-cluster",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-09T23:30:00.267496Z",
     "iopub.status.busy": "2026-10-09T23:30:00.267348Z",
     "iopub.status.idle": "2026-10-09T23:30:08.963457Z",
     "shell.execute_reply": "2026-10-09T23:30:08.962491Z"
    }
   },
   "outputs": [],
   "source": [
    "indicator = ClusterAnalysisIndicator(csv_path, coord_cols=(\"x\", \"y\", \"z\"))\n",
    "spec = SimilarityMetricLibrary().von_mises_stress()\n",
    "\n",
    "run = indicator.run(\n",
    "    method_type=\"scipy_hierarchical\",\n",
    "    spec=spec,\n",
    "    threshold=0.01,\n",
    "    method=\"average\",\n",
    "    criterion=\"distance\",\n",
    ")\n",
    "result = run[\"points\"]\n",
    "print(result[[\"x\", \"y\", \"cluster_label\"]].head())\n",
    "print(\"n clusters:\", result[\"cluster_label\"].nunique())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "r2d-viz-md",
   "metadata": {},
   "source": [
    "## Visualize field vs. clusters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "r2d-viz",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-09T23:30:08.965293Z",
     "iopub.status.busy": "2026-10-09T23:30:08.965118Z",
     "iopub.status.idle": "2026-10-09T23:30:09.249907Z",
     "shell.execute_reply": "2026-10-09T23:30:09.248277Z"
    }
   },
   "outputs": [],
   "source": [
    "def von_mises(df):\n",
    "    s = {k: df[k].to_numpy() for k in (\"sxx\", \"syy\", \"szz\", \"sxy\", \"sxz\", \"syz\")}\n",
    "    t1 = ((s[\"sxx\"] - s[\"syy\"]) ** 2 + (s[\"syy\"] - s[\"szz\"]) ** 2\n",
    "          + (s[\"szz\"] - s[\"sxx\"]) ** 2) / 2.0\n",
    "    t2 = 3.0 * (s[\"sxy\"] ** 2 + s[\"sxz\"] ** 2 + s[\"syz\"] ** 2)\n",
    "    return np.sqrt(t1 + t2)\n",
    "\n",
    "vm_grid = von_mises(result).reshape(ny, nx)\n",
    "labels = result[\"cluster_label\"].to_numpy()\n",
    "uniq = np.unique(labels)\n",
    "idx = {lab: i for i, lab in enumerate(uniq)}\n",
    "label_grid = np.vectorize(idx.get)(labels).reshape(ny, nx)\n",
    "cmap = ListedColormap(plt.get_cmap(\"tab20\")(np.linspace(0, 1, len(uniq))))\n",
    "norm = BoundaryNorm(np.arange(len(uniq) + 1) - 0.5, len(uniq))\n",
    "\n",
    "fig, ax = plt.subplots(1, 2, figsize=(9, 4))\n",
    "im0 = ax[0].imshow(vm_grid, origin=\"lower\", cmap=\"viridis\")\n",
    "ax[0].set_title(\"von-Mises stress\")\n",
    "fig.colorbar(im0, ax=ax[0], fraction=0.046, pad=0.04)\n",
    "ax[1].imshow(label_grid, origin=\"lower\", cmap=cmap, norm=norm)\n",
    "ax[1].set_title(f\"clusters (n={len(uniq)})\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "r2d-seealso",
   "metadata": {},
   "source": [
    "The high-stress patches come back as distinct clusters against the low-stress\n",
    "background.\n",
    "\n",
    "**See also**\n",
    "\n",
    "- Algorithm: :doc:`/algorithms/rare-event-identification`\n",
    "- The full pipeline (graph clustering + hierarchical merge + rare-cluster export)\n",
    "  in :doc:`rare-event-identification`\n",
    "- 3D version: :doc:`rei-example-3d`"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}