Detect modes

benchmarks/scripts/analyze.py

@@ -14,6 +14,10 @@
 from scipy.stats.mstats import hdquantiles
 
 
+precision = 0.04
+sensitivity = 0.5
+
+
 def get_bench_columns():
     return ['variant', 'elapsed', 'center', 'samples']
 
@@ -223,7 +227,7 @@ def coverage(args):
     iterate_case_dfs(args, case_coverage)
 
 
-def qrde_hd(samples, precision=0.01):
+def qrde_hd(samples):
     """
     Computes quantile-respectful density estimation based on the Harrell-Davis 
     quantile estimator. The implementation is based on the following post:
@@ -239,31 +243,105 @@ def qrde_hd(samples, precision=0.01):
     return width, height
 
 
-def hd_displot(samples, label, ax, precision=0.01):
-    width, height = qrde_hd(samples, precision)
+def extract_peaks(pdf):
+    peaks = []
+    for i in range(1, len(pdf) - 1):
+        if pdf[i - 1] < pdf[i] > pdf[i + 1]:
+            peaks.append(i)
+    return peaks
+
+
+def extract_modes(samples):
+    """
+    Extract modes from the given samples based on the lowland algorithm:
+    https://aakinshin.net/posts/lowland-multimodality-detection/ by Andrey Akinshin
+    Implementation is based on the https://github.com/AndreyAkinshin/perfolizer
+    LowlandModalityDetector class.
+    """
+    mode_ids = []
+
+    widths, heights = qrde_hd(samples)
+    peak_ids = extract_peaks(heights)
+    bin_area = 1.0 / len(heights)
+
+    x = min(samples)
+    peak_xs = []
+    peak_ys = []
+    bin_lower = [x]
+    for idx in range(len(heights)):
+        if idx in peak_ids:
+            peak_ys.append(heights[idx])
+            peak_xs.append(x + widths[idx] / 2)
+        x += widths[idx]
+        bin_lower.append(x)
+
+    def lowland_between(left_peak, right_peak):
+        left, right = left_peak, right_peak
+        min_height = min(heights[left_peak], heights[right_peak])
+        while left < right and heights[left] > min_height:
+            left += 1
+        while left < right and heights[right] > min_height:
+            right -= 1
+
+        width = bin_lower[right + 1] - bin_lower[left]
+        total_area = width * min_height
+        total_bin_area = (right - left + 1) * bin_area
+
+        if total_bin_area / total_area < sensitivity:
+            mode_ids.append(left_peak)
+            return True
+        return False
+
+    previousPeaks = [peak_ids[0]]
+    for i in range(1, len(peak_ids)):
+        currentPeak = peak_ids[i]
+        if heights[currentPeak] < heights[previousPeaks[-1]]:
+            if lowland_between(previousPeaks[-1], currentPeak):
+                previousPeaks = [currentPeak]
+        else:
+            while previousPeaks and (heights[previousPeaks[-1]] < heights[currentPeak]):
+                if lowland_between(previousPeaks[-1], currentPeak):
+                    previousPeaks = []
+                else:
+                    previousPeaks.pop()
+            previousPeaks.append(currentPeak)
+    mode_ids.append(previousPeaks[0])
+    return mode_ids
+
+
+def hd_displot(samples, label, ax):
+    widths, heights = qrde_hd(samples)
+    mode_ids = extract_modes(samples)
 
     min_sample, max_sample = min(samples), max(samples)
 
     xs = [min_sample]
     ys = [0]
 
+    peak_xs = []
+    peak_ys = []
+
     x = min(samples)
-    for idx in range(len(width)):
-        xs.append(x + width[idx] / 2)
-        ys.append(height[idx])
-        x += width[idx]
+    for idx in range(len(widths)):
+        xs.append(x + widths[idx] / 2)
+        ys.append(heights[idx])
+        if idx in mode_ids:
+            peak_ys.append(heights[idx])
+            peak_xs.append(x + widths[idx] / 2)
+        x += widths[idx]
 
     xs = xs + [max_sample]
     ys = ys + [0]
 
     ax.plot(xs, ys, label=label)
+    ax.plot(peak_xs, peak_ys, 'o', color='red')
     ax.legend()
     ax.fill_between(xs, ys, 0, alpha=0.4)
 
 
-def displot(data, ax, precision=0.01):
+def displot(data, ax):
     for variant in data:
-        hd_displot(data[variant], variant, ax, precision)
+        hd_displot(data[variant], variant, ax)
 
 
 def case_variants(pattern, algname, ct_point_name, case_df):
@@ -335,7 +413,7 @@ def extract_horizontal_space(df):
                 data[variant_name] = horizontal_df[horizontal_df['variant'] == variant_name].iloc[0]['samples']
 
             # sns.histplot(data=data, ax=ax, kde=True)
-            displot(data, ax, precision=0.03)
+            displot(data, ax)
 
             if len(horizontal_axes) > 0:
                 ax=axes[vertical_id, horizontal_id]