1515#include " Framework/ASoA.h"
1616#include " Framework/HistogramSpec.h" // only for VARIABLE_WIDTH
1717
18+ #include < array>
1819#include < cstddef>
1920#include < cstdint>
2021#include < tuple>
@@ -43,116 +44,44 @@ template <std::size_t N>
4344struct BinningPolicyBase {
4445 BinningPolicyBase (std::array<std::vector<double >, N> bins, bool ignoreOverflows = true ) : mBins (bins), mIgnoreOverflows (ignoreOverflows)
4546 {
46- static_assert (N <= 3 , " No default binning for more than 3 columns, you need to implement a binning class yourself " );
47- for (int i = 0 ; i < N; i++) {
47+ static_assert (N >= 1 , " A binning policy needs at least one axis " );
48+ for (std:: size_t i = 0 ; i < N; i++) {
4849 binning_helpers::expandConstantBinning (bins[i], mBins [i]);
4950 }
5051 }
5152
53+ // / Bins are numbered in row major order, the first axis varying fastest:
54+ // / bin = sum_d binOnAxis(d) * prod_{e < d} binsCount(e)
55+ // / For one, two and three axes this is the numbering this class has always
56+ // / produced. Folding over the axes only removes the ceiling that came from
57+ // / writing the nest out by hand; it is not a new numbering scheme.
58+ // /
59+ // / The index is used by groupTable() as an equality key and as a total order, so
60+ // / all that is required of it is injectivity. Bins that are adjacent in parameter
61+ // / space are not adjacent in index, and nothing relies on them being so.
5262 template <typename ... Ts>
5363 int getBin (std::tuple<Ts...> const & data) const
5464 {
5565 static_assert (sizeof ...(Ts) == N, " There must be the same number of binning axes and data values/columns" );
5666
57- unsigned int i = 2 , j = 2 , k = 2 ;
58- if (this ->mIgnoreOverflows ) {
59- // underflow
60- if (std::get<0 >(data) < this ->mBins [0 ][1 ]) { // mBins[0][0] is a dummy VARIABLE_WIDTH
61- return -1 ;
67+ // Fold the heterogeneous values into one array, so that the search below is a
68+ // plain loop over the axes instead of a nest whose depth has to be spelled out.
69+ // Every value was already compared against a double bin edge before, so going
70+ // through double here changes nothing.
71+ std::array<double , N> const values = std::apply (
72+ [](auto const &... value) { return std::array<double , N>{static_cast <double >(value)...}; }, data);
73+
74+ int bin = 0 ;
75+ int stride = 1 ;
76+ for (std::size_t axis = 0 ; axis < N; axis++) {
77+ unsigned int edge = 0 ;
78+ if (!findEdge (axis, values[axis], edge)) {
79+ return -1 ; // outside this axis, and outside values are dropped
6280 }
63- if constexpr (N > 1 ) {
64- if (std::get<1 >(data) < this ->mBins [1 ][1 ]) { // mBins[1][0] is a dummy VARIABLE_WIDTH
65- return -1 ;
66- }
67- }
68- if constexpr (N > 2 ) {
69- if (std::get<2 >(data) < this ->mBins [2 ][1 ]) { // mBins[2][0] is a dummy VARIABLE_WIDTH
70- return -1 ;
71- }
72- }
73- } else {
74- i = 1 ;
75- j = 1 ;
76- k = 1 ;
81+ bin += (static_cast <int >(edge) - 1 - getOverflowShift ()) * stride;
82+ stride *= getBinsCount (mBins [axis]);
7783 }
78-
79- for (; i < this ->mBins [0 ].size (); i++) {
80- if (std::get<0 >(data) < this ->mBins [0 ][i]) {
81-
82- if constexpr (N > 1 ) {
83- for (; j < this ->mBins [1 ].size (); j++) {
84- if (std::get<1 >(data) < this ->mBins [1 ][j]) {
85-
86- if constexpr (N > 2 ) {
87- for (; k < this ->mBins [2 ].size (); k++) {
88- if (std::get<2 >(data) < this ->mBins [2 ][k]) {
89- return getBinAt (i, j, k);
90- }
91- }
92- if (this ->mIgnoreOverflows ) {
93- return -1 ;
94- }
95- }
96-
97- // overflow for mBins[2] only
98- return getBinAt (i, j, k);
99- }
100- }
101-
102- if (this ->mIgnoreOverflows ) {
103- return -1 ;
104- }
105-
106- // overflow for mBins[1] only
107- if constexpr (N > 2 ) {
108- for (k = 2 ; k < this ->mBins [2 ].size (); k++) {
109- if (std::get<2 >(data) < this ->mBins [2 ][k]) {
110- return getBinAt (i, j, k);
111- }
112- }
113- }
114- }
115-
116- // overflow for mBins[2] and mBins[1]
117- return getBinAt (i, j, k);
118- }
119- }
120-
121- if (this ->mIgnoreOverflows ) {
122- // overflow
123- return -1 ;
124- }
125-
126- // overflow for mBins[0] only
127- if constexpr (N > 1 ) {
128- for (j = 2 ; j < this ->mBins [1 ].size (); j++) {
129- if (std::get<1 >(data) < this ->mBins [1 ][j]) {
130-
131- if constexpr (N > 2 ) {
132- for (k = 2 ; k < this ->mBins [2 ].size (); k++) {
133- if (std::get<2 >(data) < this ->mBins [2 ][k]) {
134- return getBinAt (i, j, k);
135- }
136- }
137- }
138-
139- // overflow for mBins[0] and mBins[2]
140- return getBinAt (i, j, k);
141- }
142- }
143- }
144-
145- // overflow for mBins[0] and mBins[1]
146- if constexpr (N > 2 ) {
147- for (k = 2 ; k < this ->mBins [2 ].size (); k++) {
148- if (std::get<2 >(data) < this ->mBins [2 ][k]) {
149- return getBinAt (i, j, k);
150- }
151- }
152- }
153-
154- // overflow for all bins
155- return getBinAt (i, j, k);
84+ return bin;
15685 }
15786
15887 // Note: Overflow / underflow bin -1 is not included
@@ -179,46 +108,53 @@ struct BinningPolicyBase {
179108 return getBinsCount (mBins [2 ]);
180109 }
181110
111+ // Note: Overflow / underflow bin -1 is not included
112+ int getBinsCountForAxis (std::size_t axis) const
113+ {
114+ return getBinsCount (mBins [axis]);
115+ }
116+
182117 // Note: Overflow / underflow bin -1 is not included
183118 int getAllBinsCount () const
184119 {
185- if constexpr (N == 1 ) {
186- return getXBinsCount ();
187- }
188- if constexpr (N == 2 ) {
189- return getXBinsCount () * getYBinsCount ();
120+ int count = 1 ;
121+ for (std::size_t axis = 0 ; axis < N; axis++) {
122+ count *= getBinsCount (mBins [axis]);
190123 }
191- if constexpr (N == 3 ) {
192- return getXBinsCount () * getYBinsCount () * getZBinsCount ();
193- }
194- return -1 ;
124+ return count;
195125 }
196126
197127 std::array<std::vector<double >, N> mBins ;
198128 bool mIgnoreOverflows ;
199129
200130 private:
201- // We substract 1 to account for VARIABLE_WIDTH in the bins vector
202- // We substract second 1 if we omit values below minima (underflow, mapped to -1)
203- // Otherwise we add 1 and we get the number of bins including those below and over the outer edges
204- int getBinAt (unsigned int iRaw, unsigned int jRaw, unsigned int kRaw ) const
131+ // / Index of the first edge strictly above the value. mBins[axis][0] is a dummy
132+ // / VARIABLE_WIDTH marker and mBins[axis][1] is the lower edge, so the first
133+ // / candidate is 2 when under- and overflows are dropped, and 1 when they are kept
134+ // / and get bins of their own. Returns false when the value falls outside the axis
135+ // / and outside values are being dropped.
136+ bool findEdge (std::size_t axis, double value, unsigned int & edge) const
205137 {
206- int shiftBinsWithoutOverflow = getOverflowShift ();
207- unsigned int i = iRaw - 1 - shiftBinsWithoutOverflow;
208- unsigned int j = jRaw - 1 - shiftBinsWithoutOverflow;
209- unsigned int k = kRaw - 1 - shiftBinsWithoutOverflow;
210- auto xBinsCount = getXBinsCount ();
211- if constexpr (N == 1 ) {
212- return i;
213- } else if constexpr (N == 2 ) {
214- return i + j * xBinsCount;
215- } else if constexpr (N == 3 ) {
216- return i + j * xBinsCount + k * xBinsCount * getYBinsCount ();
217- } else {
218- return -1 ;
138+ auto const & edges = mBins [axis];
139+ if (mIgnoreOverflows && value < edges[1 ]) {
140+ return false ; // underflow
141+ }
142+ for (unsigned int i = mIgnoreOverflows ? 2 : 1 ; i < edges.size (); i++) {
143+ if (value < edges[i]) {
144+ edge = i;
145+ return true ;
146+ }
219147 }
148+ if (mIgnoreOverflows ) {
149+ return false ; // overflow
150+ }
151+ edge = static_cast <unsigned int >(edges.size ());
152+ return true ;
220153 }
221154
155+ // We substract 1 to account for VARIABLE_WIDTH in the bins vector
156+ // We substract second 1 if we omit values below minima (underflow, mapped to -1)
157+ // Otherwise we add 1 and we get the number of bins including those below and over the outer edges
222158 int getOverflowShift () const
223159 {
224160 return mIgnoreOverflows ? 1 : -1 ;
0 commit comments