Skip to content

[Math] Fix bug and warning in recently-introduced TMath::ModeHalfSample #19578

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

[Math] Fix bug and warning in recently-introduced TMath::ModeHalfSample #19578

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
42 changes: 26 additions & 16 deletions math/mathcore/inc/TMath.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1451,23 +1451,33 @@ template <typename T> Double_t TMath::ModeHalfSample(Long64_t n, const T *a, con
if (w && std::adjacent_find(w, w + n, std::not_equal_to<>()) == w + n) // If all weights are equal, ignore those
w = nullptr;

// Sort the array and remove duplicates (if w != nullptr)
// Sort the array
std::vector<T> values;
std::vector<Double_t> weights;
values.reserve(n);
weights.reserve(n);
for (Long64_t i = 0; i < n; ++i) {
const T value = a[i];
const Double_t weight = w ? w[i] : 1;
const auto vstart = values.begin();
const auto vstop = values.end();
const auto viter = std::lower_bound(vstart, vstop, value);
const auto vidx = std::distance(vstart, viter);
if (viter != vstop && w) {
weights[vidx] += weight;
} else {
values.insert(viter, value);
weights.insert(weights.begin() + vidx, weight);
if (w == nullptr) {
values.assign(a, a + n);
weights.assign(n, 1.0);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
weights.assign(n, 1.0);

std::sort(values.begin(), values.end());
} else {
// Remove duplicates (if w != nullptr).
// The Bickel algorithm expects all values have the same weight of 1.
// The deduplication is relevant for varying weights, where we fall back
// to getting the Mode of the sample.
values.reserve(n);
weights.reserve(n);
for (Long64_t i = 0; i < n; ++i) {
const T value = a[i];
const auto vstart = values.begin();
const auto viter = std::lower_bound(vstart, values.end(), value);
const auto vidx = std::distance(vstart, viter);
if (viter != values.end() && *viter == value) {
// increase weight if value already added
weights[vidx] += w[i];
} else {
// otherwise, insert at the right position
values.insert(viter, value);
weights.insert(weights.begin() + vidx, w[i]);
}
}
}

Expand All @@ -1483,7 +1493,7 @@ template <typename T> Double_t TMath::ModeHalfSample(Long64_t n, const T *a, con
const double d1_0 = values[1] - values[0];
const double d2_1 = values[2] - values[1];
if (d2_1 < d1_0)
return TMath::Mean(2, values.data() + 1, weights.data() + 1);
return TMath::Mean(std::next(values.begin()), values.end(), std::next(weights.begin()));
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If the goal is optimization, then this should maybe also go in a separate codepath, one no weights variable, and one with.

Copy link
Collaborator

@ferdymercury ferdymercury Aug 8, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

1480: same comment, could be separate.
1505: this is a "fallback-jump" for very rare cases. In principle we were about to get the normal mode of the sample because weights were varying. But if after deduplicating and summing, weights end up being all equal, then we switch to Bickel algorithm.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please also remove "are unique" in the comment in 1512, thanks!

else if (d2_1 > d1_0)
return TMath::Mean(2, values.data(), weights.data());
else
Expand Down
Loading