{ "cells": [ { "cell_type": "markdown", "id": "2e7f55aa", "metadata": {}, "source": [ "# Rare-event identification (REI)\n", "\n", "A CPFE field has millions of numbers in it and a handful of places that matter.\n", "REI is the machinery for going from the first to the second: find the regions\n", "that are both **extreme in a field** and **spatially coherent**, and hand them\n", "back as labelled clusters you can export, count and compare.\n", "\n", "The hard part is the word *coherent*. Thresholding a field gives you the extreme\n", "points, but they come back as speckle \u2014 isolated voxels scattered through the\n", "domain, most of them noise. What you actually want is a connected region that is\n", "extreme *together*. REI gets there in two stages:\n", "\n", "1. **Over-segment** (`GraphSpatialCluster`). Build a graph on the points, weight\n", " each edge by how similar its endpoints are in the chosen field, and run Leiden\n", " at a high resolution. This produces many small, internally uniform clusters.\n", " Spatial adjacency is baked into the graph, so a cluster is connected by\n", " construction.\n", "2. **Merge in feature space** (`ClusterAnalysisIndicator`). Reduce each stage-1\n", " cluster to its mean feature vector and agglomerate hierarchically. Cutting\n", " that dendrogram at `threshold` groups clusters that are *alike*, whether or\n", " not they touch.\n", "\n", "Then `RareCriteria` picks the rare ones out of the stage-2 result. Splitting it\n", "this way is what lets the two stages ask different questions \u2014 stage 1 is \"is\n", "this one region?\", stage 2 is \"is this region like that one?\".\n", "\n", "**Needs:** pip install only (networkit Leiden, scipy hierarchical clustering,\n", "PyVista/VTK).\n", "\n", "See the algorithm page [Rare-event identification](../algorithms/rare-event-identification.rst)\n", "and the API `IdentifyRareClusters`." ] }, { "cell_type": "markdown", "id": "6b706c85", "metadata": {}, "source": [ "## Inputs\n", "\n", "A point-cloud CSV with an `id` column, `x, y, z` coordinates, and field columns.\n", "REI is not grid-locked: a regular grid uses the fast grid graph, arbitrary points\n", "fall back to the kNN graph. Two CPFE sources share the same schema:\n", "\n", "- `mesh_out/out_element_centroid_*.csv`: crisp per-element fields on the true\n", " mesh (`mesh_csv=\"sync\"`, the default), full fidelity, one row per element\n", " (kNN path). **Preferred.**\n", "- `grid_out/out_element_centroid_*.csv`: a regular grid, from\n", " `grid_transfer=\"per_step\"` or an offline `GridResampler` (smoothed; see\n", " [post-processing](post-processing.ipynb)).\n", "\n", "This notebook reads the shipped `mwe_data/synthetic_vms.csv` \u2014 a synthetic\n", "stress field with planted high-von-Mises regions \u2014 and writes everything it\n", "produces under `out/rei/`, leaving the sample data untouched." ] }, { "cell_type": "markdown", "id": "ef273326", "metadata": {}, "source": [ "## Build the metric and the stage objects\n", "\n", "A `SimilarityMetric` is the pivot of the whole pipeline: it names the feature\n", "columns and a distance function over them, and that distance is what becomes the\n", "edge weight in stage 1. Here `von_mises_stress()` takes the six stress\n", "components and defines the distance between two points as the relative\n", "difference of their von Mises equivalent stress \u2014 so two points are \"similar\"\n", "when they are carrying a similar equivalent load, regardless of the stress\n", "state's orientation.\n", "\n", "Stage 2 needs its own metric over the *reduced* table, where every feature\n", "column has picked up a `_mean` suffix during the reduction. `spec_reduced`\n", "is the same distance function pointed at those renamed columns.\n", "\n", "`WeightConfig` turns a distance into a weight. `mode=\"rbf\"` gives\n", "$w = \\exp(-(d/\\sigma)^{p})$, bounded in $(0, 1]$, which keeps the Leiden\n", "objective well conditioned; `sigma_auto` sets $\\sigma$ from a quantile of the\n", "sampled edge distances. (On near-identical *orientation* data `sigma_auto`\n", "collapses and shatters the segmentation \u2014 see\n", "[segmentation](../algorithms/segmentation.rst) \u2014 but on a smoothly varying\n", "stress field it is the right default.)" ] }, { "cell_type": "code", "execution_count": null, "id": "ac4d2644", "metadata": { "execution": { "iopub.execute_input": "2026-10-09T23:27:06.266503Z", "iopub.status.busy": "2026-10-09T23:27:06.265742Z", "iopub.status.idle": "2026-10-09T23:27:06.705666Z", "shell.execute_reply": "2026-10-09T23:27:06.704300Z" } }, "outputs": [], "source": [ "import os\n", "\n", "import numpy as np\n", "import pandas as pd\n", "\n", "from graintrace.rare_cluster_indicator import IdentifyRareClusters\n", "from graintrace.similarity_metric_library import SimilarityMetricLibrary\n", "from graintrace.user_data_class import SimilarityMetric, WeightConfig, RareCriteria\n", "from graintrace import rare_criteria_selection_library as rcs\n", "\n", "os.makedirs(\"out/rei\", exist_ok=True)\n", "filename = \"../../mwe_data/synthetic_vms.csv\"\n", "\n", "metric_lib = SimilarityMetricLibrary()\n", "spec = metric_lib.von_mises_stress()\n", "spec_reduced = SimilarityMetric(\n", " name=spec.name + \"_mean\",\n", " feature_cols=[f\"{c}_mean\" for c in spec.feature_cols],\n", " func=spec.func,\n", ")\n", "\n", "weight_cfg = WeightConfig(\n", " mode=\"rbf\",\n", " power=2.0,\n", " sigma=None,\n", " sigma_auto={\"sample_size\": 20_000, \"random_state\": 42, \"quantile\": 0.5},\n", ")\n", "\n", "print(\"stage-1 features:\", spec.feature_cols)\n", "print(\"stage-2 features:\", spec_reduced.feature_cols)\n", "print(\"points:\", len(pd.read_csv(filename)))" ] }, { "cell_type": "markdown", "id": "b3c36f60", "metadata": {}, "source": [ "## Run the two-stage clustering\n", "\n", "`graph_mode=\"grid\"` builds the graph from grid adjacency within\n", "`manhattan_radius` (an $\\ell^1$ ball: radius 2 gives 12 neighbours on this\n", "$30 \\times 30$ 2D lattice, 24 in 3D), which is far cheaper than a kNN search\n", "when the points really are on a lattice. `reduce_edges_topweights_k` then keeps\n", "only each node's 20 strongest edges.\n", "\n", "`networkit_kwargs[\"gamma\"]` is the Leiden resolution. **High** here \u2014\n", "`gamma=10.0` \u2014 because stage 1 is deliberately over-segmenting: it is better to\n", "cut a real region in half and let stage 2 glue it back together than to merge\n", "two unrelated regions that stage 2 can never separate again.\n", "\n", "`threshold=0.0005` cuts the stage-2 dendrogram with `criterion=\"distance\"`, so\n", "it is a distance in the *metric's* units, not a cluster count. `method=\"average\"`\n", "is UPGMA linkage. `dendrogram_path` writes the tree out so you can see where the\n", "cut lands rather than guessing the threshold." ] }, { "cell_type": "code", "execution_count": null, "id": "90f84471", "metadata": { "execution": { "iopub.execute_input": "2026-10-09T23:27:06.709397Z", "iopub.status.busy": "2026-10-09T23:27:06.709237Z", "iopub.status.idle": "2026-10-09T23:27:08.754555Z", "shell.execute_reply": "2026-10-09T23:27:08.752986Z" } }, "outputs": [], "source": [ "irc = IdentifyRareClusters(\n", " input_csv_path=filename, id_col=\"id\", coord_cols=(\"x\", \"y\", \"z\")\n", ")\n", "gsc, indicator = irc.make_stage_objects(graph_cluster_out=\"out/rei/reduced.csv\")\n", "\n", "bundle = irc.run_clustering(\n", " gsc=gsc,\n", " indicator=indicator,\n", " reduced_csv_path=\"out/rei/reduced.csv\",\n", " gsc_run_kwargs=dict(\n", " spec=spec,\n", " graph_mode=\"grid\",\n", " manhattan_radius=2,\n", " grid_tol=1e-6,\n", " n_jobs=1,\n", " weight_chunk_size=500_000,\n", " segmenter=\"leiden\",\n", " seed=42,\n", " weight_cfg=weight_cfg,\n", " reduce_edges_topweights_k=20,\n", " networkit_kwargs={\"gamma\": 10.0},\n", " resume_from_checkpoint=False,\n", " ),\n", " indicator_run_kwargs=dict(\n", " method_type=\"scipy_hierarchical\",\n", " spec=spec_reduced,\n", " threshold=0.0005,\n", " method=\"average\",\n", " criterion=\"distance\",\n", " dendrogram_path=\"out/rei/dendrogram.png\",\n", " ),\n", ")\n", "\n", "print(\"stage-1 clusters:\", len(np.unique(bundle[\"gsc_labels\"])))\n", "print(\"stage-2 clusters:\", len(bundle[\"indicator_clusters_df\"]))" ] }, { "cell_type": "markdown", "id": "609bba33", "metadata": {}, "source": [ "The reduction is the whole point of the two-stage design: the stage-2\n", "agglomeration runs on cluster means, not on every point, which is what makes\n", "hierarchical clustering tractable on a real field at all.\n", "\n", "Read the counts honestly. 900 points became 65 stage-1 clusters, and\n", "`threshold=0.0005` merged those into 51 \u2014 so most stage-1 clusters survived as\n", "singletons. That is a *tight* cut, and on this small synthetic field it is close\n", "to no merging at all. The dendrogram below is how you calibrate it: find the\n", "height at which the branches you consider \"the same kind of region\" join, and\n", "set `threshold` just above it. A larger threshold merges more aggressively and\n", "gives fewer, broader rare regions." ] }, { "cell_type": "code", "execution_count": null, "id": "5a283dbd", "metadata": { "execution": { "iopub.execute_input": "2026-10-09T23:27:08.758743Z", "iopub.status.busy": "2026-10-09T23:27:08.758576Z", "iopub.status.idle": "2026-10-09T23:27:08.768572Z", "shell.execute_reply": "2026-10-09T23:27:08.766597Z" } }, "outputs": [], "source": [ "from IPython.display import Image\n", "\n", "Image(\"out/rei/dendrogram.png\")" ] }, { "cell_type": "markdown", "id": "f6d9c316", "metadata": {}, "source": [ "## Select the rare clusters\n", "\n", "`bundle[\"indicator_clusters_df\"]` carries one row per stage-2 cluster with\n", "`cluster_label`, the point count `n`, and `min/max/sum/sumsq/mean` for every\n", "coordinate and every reduced feature. Note the double suffix: the stage-1\n", "reduction renamed `szz` to `szz_mean`, and the stage-2 aggregation then takes\n", "the mean of *that*, giving `szz_mean_mean`.\n", "\n", "There is **no `von_mises_stress_mean_mean` column** \u2014 the aggregation is\n", "per-component, and von Mises is a nonlinear function of the components, so the\n", "pipeline cannot produce it by averaging. Derive it in the selector instead,\n", "which is exactly what `RareCriteria` is for: it takes any callable from the\n", "cluster table to a list of labels." ] }, { "cell_type": "code", "execution_count": null, "id": "d7d07d78", "metadata": { "execution": { "iopub.execute_input": "2026-10-09T23:27:08.772273Z", "iopub.status.busy": "2026-10-09T23:27:08.772064Z", "iopub.status.idle": "2026-10-09T23:27:08.777031Z", "shell.execute_reply": "2026-10-09T23:27:08.775660Z" } }, "outputs": [], "source": [ "print([c for c in bundle[\"indicator_clusters_df\"].columns if c.endswith(\"_mean_mean\")])" ] }, { "cell_type": "code", "execution_count": null, "id": "6b379b0b", "metadata": { "execution": { "iopub.execute_input": "2026-10-09T23:27:08.780325Z", "iopub.status.busy": "2026-10-09T23:27:08.780095Z", "iopub.status.idle": "2026-10-09T23:27:08.794008Z", "shell.execute_reply": "2026-10-09T23:27:08.792078Z" } }, "outputs": [], "source": [ "def highest_von_mises(df, k=3, min_size=1):\n", " \"\"\"Rank stage-2 clusters by the von Mises stress of their mean stress tensor.\"\"\"\n", " d = df.copy()\n", " s = {c: d[f\"{c}_mean_mean\"].to_numpy() for c in spec.feature_cols}\n", " d[\"von_mises_mean\"] = np.sqrt(\n", " 0.5\n", " * (\n", " (s[\"sxx\"] - s[\"syy\"]) ** 2\n", " + (s[\"syy\"] - s[\"szz\"]) ** 2\n", " + (s[\"szz\"] - s[\"sxx\"]) ** 2\n", " )\n", " + 3.0 * (s[\"sxy\"] ** 2 + s[\"syz\"] ** 2 + s[\"sxz\"] ** 2)\n", " )\n", " return rcs.select_highest_scalar(\n", " d, k=k, required_cols=\"von_mises_mean\", min_size=min_size\n", " )\n", "\n", "\n", "rare_criteria = RareCriteria(selector=highest_von_mises)\n", "print(\"rare labels:\", highest_von_mises(bundle[\"indicator_clusters_df\"]))" ] }, { "cell_type": "markdown", "id": "62e487e1", "metadata": {}, "source": [ "## Export\n", "\n", "`run_get_rare_cluster` applies the criteria and writes a VTK with the rare\n", "clusters as labelled blocks: `background_block_id=1` for everything not\n", "selected, then `first_rare_block_id=2` upward, one id per rare cluster. Open it\n", "in ParaView and threshold on the block id to see the rare regions in place.\n", "\n", "`rare_points_csv_path` additionally writes the bare `x, y, z, rare_cluster_id`\n", "point cloud, which is the input format\n", "[REI comparison](rei-comparison.ipynb) expects \u2014 that is how you ask whether two\n", "different criteria flag the same places." ] }, { "cell_type": "code", "execution_count": null, "id": "0f105997", "metadata": { "execution": { "iopub.execute_input": "2026-10-09T23:27:08.797913Z", "iopub.status.busy": "2026-10-09T23:27:08.797381Z", "iopub.status.idle": "2026-10-09T23:27:08.856508Z", "shell.execute_reply": "2026-10-09T23:27:08.855278Z" } }, "outputs": [], "source": [ "out = irc.run_get_rare_cluster(\n", " bundle=bundle,\n", " criteria=rare_criteria,\n", " output_vtk_path=\"out/rei/rare_clusters.vtk\",\n", " export_control=\"auto\",\n", " background_block_id=1,\n", " first_rare_block_id=2,\n", " also_write_final_label=True,\n", " rare_reduced_stats_csv_path=\"out/rei/rare_cluster_stats.csv\",\n", " rare_points_csv_path=\"out/rei/rare_points.csv\",\n", ")\n", "print(\"VTK exported:\", out[\"output_vtk_path\"])\n", "\n", "stats = pd.read_csv(\"out/rei/rare_cluster_stats.csv\")\n", "print(stats[[\"cluster_label\", \"n\", \"x_mean\", \"y_mean\", \"z_mean\"]].to_string(index=False))" ] }, { "cell_type": "markdown", "id": "f9b98382", "metadata": {}, "source": [ "## Notes\n", "\n", "- **Checkpointing.** Three restart levels exist for long runs: the bundle\n", " pickle, the reduced CSV plus the stage-1 labels array, and the graph edges\n", " themselves (`checkpoint_base_path` + `resume_from_checkpoint=True`). The graph\n", " checkpoint deliberately does **not** compare `segmenter` or the NetworKit\n", " options, so one checkpoint is reusable across a whole gamma sweep. See\n", " [configuration](../configuration.rst).\n", "- **No multiprocessing.** The pipeline is numba-threaded and single-process, so\n", " no `if __name__ == \"__main__\":` guard is required \u2014 unlike NF reconstruction.\n", "- **`sigma_auto` is field-dependent.** Fine on a smooth stress field, actively\n", " harmful on orientation data where intra-grain edges are near-identical.\n", "\n", "**See also**\n", "\n", "- Algorithm: [Rare-event identification](../algorithms/rare-event-identification.rst)\n", "- Smaller worked examples: [2D](rei-example-2d.ipynb), [3D](rei-example-3d.ipynb)\n", "- Compare two REI outputs: [REI comparison](rei-comparison.ipynb)\n", "- Full script: [examples/demonstrate_rei_pipeline.py](https://github.com/applied-material-modeling/graintrace/blob/main/examples/demonstrate_rei_pipeline.py)" ] } ], "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 }