From 0de70dde7e23a6718c2ff35284db470f338c303a Mon Sep 17 00:00:00 2001
From: Filip Leonarski <filip.leonarski@psi.ch>
Date: Fri, 12 Dec 2025 14:28:06 +0100
Subject: [PATCH] Histogram: Add Percentile calculation

---
 common/Histogram.h      | 59 +++++++++++++++++++++++++++++++++++------
 tests/HistogramTest.cpp | 17 ++++++++++++
 2 files changed, 68 insertions(+), 8 deletions(-)

diff --git a/common/Histogram.h b/common/Histogram.h
index 4831aec1..dc33d0d8 100644
--- a/common/Histogram.h
+++ b/common/Histogram.h
@@ -12,13 +12,16 @@
 #include "MultiLinePlot.h"
 #include "../common/JFJochException.h"
 
-template <class T>
+template<class T>
 class SetAverage {
     std::vector<T> sum;
     std::vector<uint64_t> count;
-    std::mutex m;
+    mutable std::mutex m;
+
 public:
-    explicit SetAverage(size_t bins) : sum(bins), count(bins) {}
+    explicit SetAverage(size_t bins) : sum(bins), count(bins) {
+    }
+
     void Add(size_t bin, T val) {
         std::unique_lock ul(m);
 
@@ -28,7 +31,7 @@ public:
         }
     }
 
-    MultiLinePlot GetPlot() {
+    MultiLinePlot GetPlot() const {
         std::unique_lock ul(m);
 
         MultiLinePlotStruct plot;
@@ -49,16 +52,19 @@ public:
 
 class FloatHistogram {
     std::vector<uint64_t> count;
-    std::mutex m;
+    mutable std::mutex m;
     float min;
     float max;
     float spacing;
+
 public:
     explicit FloatHistogram(size_t bins, float in_min, float in_max) : count(bins), min(in_min), max(in_max) {
         if (bins == 0)
-            throw JFJochException(JFJochExceptionCategory::InputParameterBelowMin, "FloatHistogram neds to have non-zero bins");
+            throw JFJochException(JFJochExceptionCategory::InputParameterBelowMin,
+                                  "FloatHistogram neds to have non-zero bins");
         spacing = (max - min) / bins;
     }
+
     void Add(float val) {
         std::unique_lock ul(m);
 
@@ -68,11 +74,11 @@ public:
             count[bin] += 1;
     }
 
-    MultiLinePlot GetPlot() {
+    MultiLinePlot GetPlot() const {
         std::unique_lock ul(m);
 
         MultiLinePlotStruct plot;
-       plot.x.resize(count.size());
+        plot.x.resize(count.size());
         plot.y.resize(count.size());
         for (int i = 0; i < count.size(); i++) {
             plot.x[i] = min + (i + 0.5f) * spacing;
@@ -82,6 +88,43 @@ public:
         ret.AddPlot(plot);
         return ret;
     }
+
+    uint64_t TotalCount() const {
+        std::unique_lock ul(m);
+        uint64_t total = 0;
+        for (auto c: count) total += c;
+        return total;
+    }
+
+    // Returns the value x such that approximately `percent`% of samples are <= x.
+    // - percent must be in [0, 100]
+    // - returns std::nullopt if histogram is empty
+    std::optional<float> Percentile(float percent) const {
+        if (!std::isfinite(percent) || percent < 0.0f || percent > 100.0f) {
+            throw JFJochException(JFJochExceptionCategory::InputParameterBelowMin,
+                                  "FloatHistogram Percentile expects percent in [0, 100]");
+        }
+
+        std::unique_lock ul(m);
+
+        uint64_t total = 0;
+        for (auto c: count) total += c;
+        if (total == 0) return std::nullopt;
+
+        // Target rank in [0, total-1]
+        const double q = static_cast<double>(percent) / 100.0;
+        const auto target = static_cast<uint64_t>(std::floor(q * static_cast<double>(total - 1)));
+
+        uint64_t cumulative = 0;
+        for (size_t i = 0; i < count.size(); i++) {
+            cumulative += count[i];
+            if (target < cumulative && count[i] > 0)
+                    return min + static_cast<float>(i) * spacing + 0.5f * spacing;
+        }
+
+        // If due to rounding we didn't return inside the loop, clamp to the last bin's upper edge.
+        return max;
+    }
 };
 
 #endif //JUNGFRAUJOCH_HISTOGRAM_H
diff --git a/tests/HistogramTest.cpp b/tests/HistogramTest.cpp
index 22ab5c26..b808a1a4 100644
--- a/tests/HistogramTest.cpp
+++ b/tests/HistogramTest.cpp
@@ -53,3 +53,20 @@ TEST_CASE("FloatHistogram") {
     CHECK(p.GetPlots()[0].y[50] == 1.0);
     CHECK(p.GetPlots()[0].y[99] == 0.0);
 }
+
+TEST_CASE("FloatHistogram_Percentile") {
+    FloatHistogram h(111, 89.5, 200.5);
+    for (int i = 0; i < 8; i++)
+        h.Add(100.0f);
+    h.Add(150.0f);
+    h.Add(200.0f);
+
+    CHECK(h.Percentile(0.0f).value_or(-1) == Catch::Approx(100.0f));
+    CHECK(h.Percentile(30.0f).value_or(-1) == Catch::Approx(100.0f));
+    CHECK(h.Percentile(50.0f).value_or(-1) == Catch::Approx(100.0f));
+    CHECK(h.Percentile(80.0f).value_or(-1) == Catch::Approx(100.0f));
+    CHECK(h.Percentile(90.0f).value_or(-1) == Catch::Approx(150.0f));
+    CHECK(h.Percentile(100.0f).value_or(-1) == Catch::Approx(200.0f));
+
+    CHECK(h.TotalCount() == 10);
+}