@@ -4975,71 +4975,122 @@ struct VM {
49754975 std::vector<u64 > s = samples_in;
49764976 std::sort (s.begin (), s.end ()); // ascending
49774977
4978- // trivial small -sample handling
4978+ // tiny -sample short-circuits
49794979 if (N <= 4 ) return s.front ();
49804980
4981- // Compute gaps between consecutive sorted samples
4982- std::vector<u64 > gaps;
4983- gaps.reserve (N - 1 );
4984- for (size_t i = 1 ; i < N; ++i) gaps.push_back (s[i] - s[i - 1 ]);
4985-
4986- // median gap via nth_element
4987- std::vector<u64 > gaps_copy = gaps;
4988- const size_t mid_idx = gaps_copy.size () / 2 ;
4989- std::nth_element (gaps_copy.begin (), gaps_copy.begin () + static_cast <std::ptrdiff_t >(mid_idx), gaps_copy.end ());
4990- const u64 median_gap = gaps_copy[mid_idx];
4991-
4992- // heuristics / parameters
4993- constexpr double GAP_FACTOR = 5.0 ; // require gap >= GAP_FACTOR * median_gap
4994- constexpr u64 GAP_ABS_MIN = 50 ; // or an absolute minimum gap
4995- constexpr double LOW_PERCENTILE = 0.10 ; // fallback if no gap found
4996- constexpr double TRIM_RATIO = 0.10 ; // trimmed mean ratio inside cluster
4997- constexpr double MIN_CLEAN_FRACTION = 0.05 ; // require at least this fraction to accept cluster
4998-
4999- const u64 gap_threshold = static_cast <u64 >(std::max<double >(static_cast <double >(GAP_ABS_MIN ), std::ceil (GAP_FACTOR * static_cast <double >(median_gap))));
5000-
5001- // find first "large" gap and split index is i+1 (samples[0..i] is low cluster)
5002- size_t split_index = 0 ;
5003- for (size_t i = 0 ; i < gaps.size (); ++i) {
5004- if (gaps[i] >= gap_threshold) { split_index = i + 1 ; break ; }
5005- }
5006-
5007- // fallback to low-percentile if no clear gap
5008- if (split_index == 0 ) {
5009- split_index = static_cast <size_t >(std::max<size_t >(1 , static_cast <size_t >(std::floor (static_cast <double >(N) * LOW_PERCENTILE ))));
5010- }
5011-
5012- if (split_index > N) split_index = N;
5013- size_t cluster_size = split_index;
5014-
5015- // if cluster is too small relative to N, use percentile fallback
5016- if (static_cast <double >(cluster_size) / static_cast <double >(N) < MIN_CLEAN_FRACTION ) {
5017- cluster_size = static_cast <size_t >(std::max<size_t >(1 , static_cast <size_t >(std::floor (static_cast <double >(N) * LOW_PERCENTILE ))));
5018- if (cluster_size > N) cluster_size = N;
5019- }
4981+ // median (and works for sorted input)
4982+ auto median_of_sorted = [](const std::vector<u64 >& v, size_t lo, size_t hi) -> u64 {
4983+ // this is the median of v[lo..hi-1], requires 0 <= lo < hi
4984+ const size_t len = hi - lo;
4985+ if (len == 0 ) return 0 ;
4986+ const size_t mid = lo + (len / 2 );
4987+ if (len & 1 ) return v[mid];
4988+ return (v[mid - 1 ] + v[mid]) / 2 ;
4989+ };
50204990
5021- // compute robust estimate for cluster s[0 ... cluster_size-1]
5022- u64 result = 0 ;
5023- if (cluster_size >= 10 ) {
5024- const size_t trim = static_cast <size_t >(std::floor (static_cast <double >(cluster_size) * TRIM_RATIO ));
5025- const size_t lo = trim;
5026- const size_t hi = cluster_size - trim; // exclusive
5027- if (hi <= lo) {
5028- // degenerate which is cluster median
5029- const size_t mid = cluster_size / 2 ;
5030- result = (cluster_size % 2 ) ? s[mid] : ((s[mid - 1 ] + s[mid]) / 2 );
5031- }
5032- else {
5033- unsigned long long sum = 0 ;
5034- for (size_t i = lo; i < hi; ++i) sum += s[i];
5035- result = static_cast <u64 >(static_cast <double >(sum) / static_cast <double >(hi - lo) + 0.5 );
5036- }
5037- }
5038- else {
5039- // small cluster which is median of cluster
5040- const size_t mid = cluster_size / 2 ;
5041- result = (cluster_size % 2 ) ? s[mid] : ((s[mid - 1 ] + s[mid]) / 2 );
5042- }
4991+ // the robust center: median M and MAD -> approximate sigma
4992+ const u64 M = median_of_sorted (s, 0 , s.size ());
4993+ std::vector<u64 > absdev;
4994+ absdev.reserve (N);
4995+ for (size_t i = 0 ; i < N; ++i) {
4996+ const u64 d = (s[i] > M) ? (s[i] - M) : (M - s[i]);
4997+ absdev.push_back (d);
4998+ }
4999+ std::sort (absdev.begin (), absdev.end ());
5000+ const u64 MAD = median_of_sorted (absdev, 0 , absdev.size ());
5001+ // convert MAD to an approximate standard-deviation-like measure
5002+ const long double kMADtoSigma = 1 .4826L ; // consistent for normal approx
5003+ const long double sigma = (MAD == 0 ) ? 1 .0L : (static_cast <long double >(MAD ) * kMADtoSigma );
5004+
5005+ // find the densest small-valued cluster by sliding a fixed-count window
5006+ // this locates the most concentrated group of samples (likely it would be the true VMEXIT cluster)
5007+ // const size_t frac_win = (N * 8 + 99) / 100; // ceil(N * 0.08)
5008+ // const size_t win = std::min(N, std::max(MIN_WIN, frac_win));
5009+ const size_t MIN_WIN = 10 ;
5010+ const size_t win = std::min (
5011+ N,
5012+ std::max (
5013+ MIN_WIN ,
5014+ static_cast <size_t >(std::ceil (static_cast <double >(N) * 0.08 ))
5015+ )
5016+ );
5017+ size_t best_i = 0 ;
5018+ u64 best_span = (s.back () - s.front ()) + 1 ; // large initial
5019+ for (size_t i = 0 ; i + win <= N; ++i) {
5020+ const u64 span = s[i + win - 1 ] - s[i];
5021+ if (span < best_span) {
5022+ best_span = span;
5023+ best_i = i;
5024+ }
5025+ }
5026+
5027+ // expand the initial window greedily while staying "tight"
5028+ // allow expansion while adding samples does not more than multiply the span by EXPAND_FACTOR
5029+ constexpr long double EXPAND_FACTOR = 1 .5L ;
5030+ size_t cluster_lo = best_i;
5031+ size_t cluster_hi = best_i + win; // exclusive
5032+ // expand left
5033+ while (cluster_lo > 0 ) {
5034+ const u64 new_span = s[cluster_hi - 1 ] - s[cluster_lo - 1 ];
5035+ if (static_cast <long double >(new_span) <= EXPAND_FACTOR * static_cast <long double >(best_span) ||
5036+ (s[cluster_hi - 1 ] <= (s[cluster_lo - 1 ] + static_cast <u64 >(std::ceil (3 .0L * sigma))))) {
5037+ --cluster_lo;
5038+ best_span = std::min (best_span, new_span);
5039+ }
5040+ else break ;
5041+ }
5042+ // expand right
5043+ while (cluster_hi < N) {
5044+ const u64 new_span = s[cluster_hi] - s[cluster_lo];
5045+ if (static_cast <long double >(new_span) <= EXPAND_FACTOR * static_cast <long double >(best_span) ||
5046+ (s[cluster_hi] <= (s[cluster_lo] + static_cast <u64 >(std::ceil (3 .0L * sigma))))) {
5047+ ++cluster_hi;
5048+ best_span = std::min (best_span, new_span);
5049+ }
5050+ else break ;
5051+ }
5052+
5053+ const size_t cluster_size = (cluster_hi > cluster_lo) ? (cluster_hi - cluster_lo) : 0 ;
5054+
5055+ // cluster must be reasonably dense and cover a non-negligible portion of samples, so this is pure sanity checks
5056+ const double fraction_in_cluster = static_cast <double >(cluster_size) / static_cast <double >(N);
5057+ const size_t MIN_CLUSTER = std::min (static_cast <size_t >(std::max<int >(5 , static_cast <int >(N / 50 ))), N); // at least 2% or 5 elements
5058+ if (cluster_size < MIN_CLUSTER || fraction_in_cluster < 0.02 ) {
5059+ // low-percentile (10th) trimmed median
5060+ const size_t fallback_count = std::max<size_t >(1 , static_cast <size_t >(std::floor (static_cast <double >(N) * 0.10 )));
5061+ // median of lowest fallback_count elements (if fallback_count==1 that's smallest)
5062+ if (fallback_count == 1 ) return s.front ();
5063+ const size_t mid = fallback_count / 2 ;
5064+ if (fallback_count & 1 ) return s[mid];
5065+ return (s[mid - 1 ] + s[mid]) / 2 ;
5066+ }
5067+
5068+ // now we try to get a robust estimate inside the cluster, trimmed mean (10% trim) centered on cluster
5069+ const size_t trim_count = static_cast <size_t >(std::floor (static_cast <double >(cluster_size) * 0.10 ));
5070+ size_t lo = cluster_lo + trim_count;
5071+ size_t hi = cluster_hi - trim_count; // exclusive
5072+ if (hi <= lo) {
5073+ // degenerate -> median of cluster
5074+ return median_of_sorted (s, cluster_lo, cluster_hi);
5075+ }
5076+
5077+ // sum with long double to avoid overflow and better rounding
5078+ long double sum = 0 .0L ;
5079+ for (size_t i = lo; i < hi; ++i) sum += static_cast <long double >(s[i]);
5080+ const long double avg = sum / static_cast <long double >(hi - lo);
5081+ u64 result = static_cast <u64 >(std::llround (avg));
5082+
5083+ // final sanity adjustments:
5084+ // if the computed result is suspiciously far from the global median (e.g., > +6*sigma)
5085+ // clamp toward the median to avoid choosing a high noisy cluster by mistake
5086+ const long double diff_from_med = static_cast <long double >(result) - static_cast <long double >(M);
5087+ if (diff_from_med > 0 && diff_from_med > (6 .0L * sigma)) {
5088+ // clamp to median + 4*sigma (conservative)
5089+ result = static_cast <u64 >(std::llround (static_cast <long double >(M) + 4 .0L * sigma));
5090+ }
5091+
5092+ // Also, if result is zero (shouldn't be) or extremely small, return a smallest observed sample
5093+ if (result == 0 ) result = s.front ();
50435094
50445095 return result;
50455096 };
0 commit comments