From 02ab470c025e0093ad9a666b222aba756cf5b9eb Mon Sep 17 00:00:00 2001 From: Rudolf Weeber Date: Mon, 15 Jun 2026 10:46:52 +0200 Subject: [PATCH] cluster_analysis: raise on fractal dimension of degenerate cluster (bug-sweep #45) Cluster::fractal_dimension fits a line through (log diameter, log particle count) points accumulated over radial shells. A 2-particle cluster (the normal minimum of pair-based clustering) yields at most one such point: with dr=0 the loop appends a single point, and with dr large enough that no second shell forms, no points are appended at all. Passing n<2 points to gsl_fit_linear, together with the residual division by log_diameters.size() (which is zero when the vector is empty), produced a silent [NaN, NaN] result with no diagnostic. The fractal dimension is mathematically undefined for fewer than two radial shells, so fabricating a value is wrong. Prevent the degenerate input by raising a clear std::runtime_error (surfacing to Python as RuntimeError via the script-interface dispatch, matching the existing sanity_checks() and no-GSL conventions) before calling gsl_fit_linear. Adds testsuite/python/cluster_analysis.py coverage that a 2-particle cluster raises RuntimeError for both dr=0 and a large dr. Co-Authored-By: Claude Opus 4.8 --- src/core/cluster_analysis/Cluster.cpp | 6 ++++++ testsuite/python/cluster_analysis.py | 21 +++++++++++++++++++++ 2 files changed, 27 insertions(+) diff --git a/src/core/cluster_analysis/Cluster.cpp b/src/core/cluster_analysis/Cluster.cpp index f718f71dce..68e00c7f89 100644 --- a/src/core/cluster_analysis/Cluster.cpp +++ b/src/core/cluster_analysis/Cluster.cpp @@ -177,6 +177,12 @@ Cluster::fractal_dimension([[maybe_unused]] double dr) { log_pcounts.push_back(log(static_cast(subcluster_ids.size()))); log_diameters.push_back(log(current_rg * 2.0)); } + if (log_diameters.size() < 2) { + throw std::runtime_error( + "The fractal dimension of a cluster with fewer than two distinct " + "radial shells (e.g. a 2-particle cluster) is undefined. Provide a " + "larger cluster or a smaller dr."); + } double c0, c1, cov00, cov01, cov11, sumsq; gsl_fit_linear(log_diameters.data(), 1, log_pcounts.data(), 1, log_pcounts.size(), &c0, &c1, &cov00, &cov01, &cov11, &sumsq); diff --git a/testsuite/python/cluster_analysis.py b/testsuite/python/cluster_analysis.py index 207fdb96a0..6b0a6689bf 100644 --- a/testsuite/python/cluster_analysis.py +++ b/testsuite/python/cluster_analysis.py @@ -146,6 +146,27 @@ def test_single_cluster_analysis_lees_edwards(self): with self.assertRaises(RuntimeError): cluster.radius_of_gyration() + @ut.skipUnless(espressomd.has_features("GSL"), + "Fractal dimension calculation requires GSL") + def test_fractal_dimension_small_cluster(self): + # A 2-particle cluster yields fewer than two radial shells, so the + # fractal dimension is mathematically undefined. The previous behaviour + # silently returned NaN (GSL linear fit with n<2 points); the call must + # raise a clear RuntimeError instead. + self.set_two_clusters() + dc = espressomd.pair_criteria.DistanceCriterion(cut_off=0.12) + self.cs.set_params(pair_criterion=dc) + self.cs.run_for_all_pairs() + + clusters = [c for _, c in self.cs.clusters] + small_cluster, = [c for c in clusters if c.size() == 2] + + # dr=0. gives a single fit point, dr large enough gives zero fit points; + # both are degenerate and must raise rather than return NaN. + for dr in (0., 10.): + with self.assertRaises(RuntimeError): + small_cluster.fractal_dimension(dr=dr) + def test_single_cluster_analysis(self): # Place particles on a line (crossing periodic boundaries) for x in np.arange(-0.2, 0.21, 0.01):