Rare-event identification#

Overview#

Rare-event identification (REI) locates spatially coherent regions of a CPFE field that are extreme in some quantity (e.g. high Nye-tensor norm, high von Mises stress). It is driven by IdentifyRareClusters, which runs a two-stage clustering: an over-segmenting graph pass (GraphSpatialCluster) followed by a hierarchical merge in feature space (ClusterAnalysisIndicator), then selects the rare merged clusters and exports them to VTK.

Why two stages#

A single clustering cannot answer the question REI asks. Partitioning the field purely in feature space finds the extreme values but scatters them across the domain, so the result is a histogram tail rather than a region. Partitioning purely in space finds compact blobs that need not be extreme in anything. REI therefore runs one pass of each, in that order: an over-segmenting graph pass constrained by the spatial graph, which guarantees every stage-1 cluster is a connected neighbourhood, then a feature-space merge over the reduced clusters, which is what makes the second pass affordable – \(10^6\) field points become \(10^3\)–\(10^4\) cluster means before any \(O(n^2)\) work happens.

The cost asymmetry is the design constraint. Stage 2 evaluates a pairwise distance matrix with a Python-callable metric, so its work grows as the square of the number of stage-1 clusters; \(\gamma\) in stage 1 is therefore as much a budget knob as a physical one.

Method#

Points are compared through a similarity metric on their field features (SimilarityMetricLibrary). Every metric is a distance: smaller means more alike, and the weight kernels of Grain segmentation invert it. Built-in metrics are

Metric	\(d(\mathbf{u}, \mathbf{v})\)	Feature columns
`nye_tensor_norm`	\(\lVert \mathbf{A}_i - \mathbf{A}_j \rVert_F\)	the nine `nye_tensor_ij`
`von_mises_stress`	relative difference, Eq. 29	six stress components
`misorientation`	disorientation angle, Eq. 17	three orientation components
`abs_scalar_diff`	\(\lvert u_1 - v_1 \rvert\)	exactly one column

The tensor metric is the Frobenius norm of the difference, not the difference of the norms, so it separates two points whose Nye tensors have equal magnitude but different structure. The von Mises metric is relative, not absolute:

(29)#\[ d = \frac{\lvert \sigma^{vM}_i - \sigma^{vM}_j\rvert} {\lvert \sigma^{vM}_i\rvert + \lvert \sigma^{vM}_j\rvert + \varepsilon}, \qquad \sigma^{vM} = \sqrt{\tfrac{1}{2}\!\left[(\sigma_{xx}-\sigma_{yy})^2 + (\sigma_{yy}-\sigma_{zz})^2 + (\sigma_{zz}-\sigma_{xx})^2\right] + 3(\sigma_{xy}^2+\sigma_{yz}^2+\sigma_{xz}^2)},\]

with \(\varepsilon = 10^{-8}\) fixed in the implementation. Normalising by the sum of the two magnitudes makes \(d\) dimensionless and bounded on \([0, 1]\), so a stage-2 threshold carries over between load steps and between materials. The cost is that it is scale-free: a 1 MPa step beside a 2 MPa neighbour is as distant as 500 beside 1000. If the absolute level is what matters, drive the run from abs_scalar_diff on a precomputed von Mises column instead.

Stage 1 (over-segment). A spatial graph over the field points (grid connectivity at a Manhattan radius, or mutual kNN) carries the metric distance on each edge, mapped to a weight by one of the kernels in Grain segmentation, and is partitioned by Leiden into many small, spatially compact clusters. Each cluster \(C\) is reduced to

(30)#\[n_C = |C|, \qquad \bar{x}_C = \frac{1}{n_C}\sum_{i \in C} x_i, \qquad \overline{f}_C = \frac{1}{n_C}\sum_{i \in C} f_i,\]

one <col>_mean per feature column. When the feature columns contain a complete \(3 \times 3\) set <prefix>_11 … <prefix>_33, an extra <prefix>_norm_mean column carries \(\tfrac{1}{n_C}\sum_{i \in C} \lVert \mathbf{A}_i \rVert_F\) – the mean of the per-point norms, not the norm of the mean tensor. The distinction matters: the mean tensor of a cluster straddling a dislocation wall partially cancels, while the mean norm does not, so rarity survives the reduction.

Stage 2 (merge). The reduced per-cluster means are agglomerated by SciPy hierarchical linkage. With the default method="average" (UPGMA) the distance between two merged clusters is the mean over all cross pairs,

(31)#\[D(\mathcal{A}, \mathcal{B}) = \frac{1}{|\mathcal{A}|\,|\mathcal{B}|} \sum_{a \in \mathcal{A}} \sum_{b \in \mathcal{B}} d(a, b),\]

which is the right default here because it is insensitive to a single outlying cluster in a way single (chaining) and complete (fragmenting) are not. The dendrogram is cut at criterion="distance" and height threshold, so threshold is in the units of the metric – dimensionless for the relative von Mises metric, degrees or radians for misorientation, and in the field’s own units for a Frobenius norm. There is no normalisation step, so a threshold transplanted between metrics is meaningless.

Quality is reported as the cophenetic correlation: the Pearson correlation between the original pairwise distances and the dendrogram heights at which each pair first merges,

(32)#\[c = \frac{\sum_{i<j} \bigl(d_{ij} - \bar{d}\bigr)\bigl(t_{ij} - \bar{t}\bigr)} {\sqrt{\sum_{i<j} \bigl(d_{ij} - \bar{d}\bigr)^2 \sum_{i<j} \bigl(t_{ij} - \bar{t}\bigr)^2}},\]

with \(t_{ij}\) the cophenetic (ultrametric) distance. A value near 1 means the tree is a faithful summary of the distance matrix and a height cut is meaningful; a low value means the data does not have a hierarchical structure and the threshold is close to arbitrary. Treat it as a go/no-go on the cut, not as a score to optimise.

Note

Stage 2 evaluates the metric on all pairs three times per run – once for the reported cophenetic_correlation, once inside linkage, and once more inside fclusterdata, which recomputes an identical tree internally to produce the labels – and then runs a 1-D MDS embedding on the dense \(n \times n\) ultrametric. All four steps are \(O(n^2)\) in the number of stage-1 clusters, with a Python call per pair. Keeping stage 1 to \(\mathcal{O}(10^3)\) clusters is what keeps stage 2 interactive.

Stage 3 (select). A RareCriteria picks which merged clusters are rare. With an explicit selector the choice is entirely the caller’s – e.g. the top-\(k\) by mean Nye-norm via graintrace.rare_criteria_selection_library.select_highest_scalar().

Warning

The default, used when selector is None, does not look at the field at all. It keeps the clusters whose point count \(n\) falls at or below the size_quantile of the surviving size distribution, after dropping anything smaller than min_size:

(33)#\[\mathcal{R} = \bigl\{\, C : \mathrm{min\_size} \le n_C \le Q_{q}(n) \,\bigr\},\]

truncated to the max_rare smallest. Its notion of rare is small, not extreme. That is a reasonable proxy only when the merge has already isolated the extreme material into a few small clusters; it is wrong whenever a rare region is large. Pass a selector for anything field-driven.

Selected clusters receive distinct block ids (background id first, then one per rare cluster, smallest first) and are written to a STRUCTURED_GRID or POLYDATA VTK, optionally with a rare point-cloud CSV for downstream comparison with REI comparison.

Algorithm#

Load the field CSV (a mesh_out/ true-mesh or grid_out/ regular-grid file) with id and coordinate columns.
Stage 1: build the spatial graph, compute per-edge metric distances, RBF-weight them, optionally prune to top-\(k\) edges per node, and Leiden-partition into fine clusters; write the reduced per-cluster CSV and per-point stage-1 labels.
Stage 2: run hierarchical linkage on the reduced cluster feature means and cut at threshold (criterion distance) to obtain merged super-labels; optionally save a dendrogram.
Map each stage-1 label to its merged super-label to get a per-point final_label.
Stage 3: apply the rarity criteria to the merged clusters, assign block ids (smallest rare clusters first), and mark rare points.
Export the classified field to VTK (grid vs. points auto-detected), plus optional per-rare-cluster statistics and an (x, y, z, rare_cluster_id) point cloud.

Parameters that matter#

spec (the SimilarityMetric): which field quantity defines rarity (Nye norm, von Mises, misorientation); stage 2 uses the *_mean reduced version.
Stage-1 graph knobs: graph_mode, manhattan_radius, weight_cfg, reduce_edges_topweights_k, and networkit_kwargs={"gamma": ...} (lower \(\gamma\) = fewer, larger clusters).
Stage-2 merge knobs: threshold (linkage cut distance), method (linkage, e.g. average), criterion (distance).
RareCriteria: selector (custom top-\(k\)) or the size_quantile / min_size / max_rare size-based default.
export_control and block-id options control the VTK output form.

See Configuration reference for the metric, weight, and criteria dataclasses.

References#

Leiden community detection (stage 1): Traag, Waltman, van Eck (2019), https://doi.org/10.1038/s41598-019-41695-z.
Agglomerative hierarchical clustering with average linkage and cophenetic distance (stage 2), as implemented in scipy.cluster.hierarchy.
The RBF edge weighting is a standard Gaussian affinity of the feature-space distance.