ATLAS Offline Software
KLGaussianMixtureReduction.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 
7 //
8 #include "CxxUtils/restrict.h"
9 #include "CxxUtils/vec.h"
10 //
11 #include <cmath>
12 #include <memory>
13 #include <limits>
14 #include <numeric>
15 #include <stdexcept>
16 #include <vector>
17 
25 namespace KLReductionFMV {
26 #ifdef __clang__ // trapping-math is incompatible with multiversioning with clang 20.1.5
27 #pragma float_control(except, off)
28 #endif
29 //clang FMV needs a namespace :/
30 #if HAVE_FUNCTION_MULTIVERSIONING
31 [[gnu::target("avx2")]]
32 int vIdxOfMin(const float* distancesIn, int n) {
33  return vAlgs::vIdxOfMin<256>(distancesIn, n);
34 }
35 [[gnu::target("default")]]
36 #endif
37 int vIdxOfMin(const float* distancesIn, int n) {
38  return vAlgs::vIdxOfMin<128>(distancesIn, n);
39 }
40 } // namespace KLReductionFMV
41 
42 namespace {
43 //internal implementation methods
44 
45 //We want to be using up to a 256 ISA, these cover also a narrower one
46 constexpr size_t STRIDEForKL = vAlgs::strideOfNumSIMDVec<256,float>(4);
47 using namespace GSFUtils;
48 
61 inline float
62 symmetricKL(const Component1D& ATH_RESTRICT componentI,
63  const Component1D& ATH_RESTRICT componentJ)
64 {
65  const double meanDifference = componentI.mean - componentJ.mean;
66  const double inverCovSum = componentI.invCov + componentJ.invCov;
67  const double term1 = componentI.invCov * componentJ.cov;
68  const double term2 = componentJ.invCov * componentI.cov;
69  const double term3 = meanDifference * inverCovSum * meanDifference;
70  return term1 + term2 + term3;
71 }
79 inline void
82 {
83 
84  const double sumWeight = updated.weight + removed.weight;
85 
86  const double invSumWeight = 1. / sumWeight;
87  const double weightI_IJ = updated.weight * invSumWeight;
88  const double weightJ_IJ = removed.weight * invSumWeight;
89  const double meanDiff = (updated.mean - removed.mean);
90 
91  const double sumMean = weightI_IJ * updated.mean + weightJ_IJ * removed.mean;
92 
93  const double sumVariance = weightI_IJ * updated.cov +
94  weightJ_IJ * removed.cov +
95  weightI_IJ * weightJ_IJ * meanDiff * meanDiff;
96  updated.mean = sumMean;
97  updated.cov = sumVariance;
98  updated.invCov = 1. / sumVariance;
99  updated.weight = sumWeight;
100 }
101 
102 inline int
103 numDistances(const int n, float* distancesIn)
104 {
105  const int npadded = vAlgs::numPadded<STRIDEForKL>(n);
106  // Make sure the padded elements are set to max
107  std::fill( distancesIn + n, distancesIn + npadded, std::numeric_limits<float>::max());
108  return npadded;
109 }
110 
111 /*
112  * This is a O(N^3) time N^2 space algorithm.
113  * It can be seen as tailored implementation for
114  * our problem of the basic algorithm
115  * for Hierarchical Clustering.
116  *
117  * We rely on the fast findMinimunIndex above
118  * and a triangular array representation to
119  * reduce the actual time for our problem
120  * (max N=72)
121  *
122  * We opt for fixed size arrays of max N(=72) elements,
123  * but dynamic arrays of ~ N*(N-1)/2 elements.
124  *
125  * Existing alternatives in the literature:
126  * - 1. O(N^3) worst-case time , O(n^2) best
127  * - 2. O(N^2 log(N)) worst-time
128  * See :
129  * "Modern hierarchical, agglomerative clustering algorithms"
130  * https://arxiv.org/abs/1109.2378
131  * or
132  * "Efficient algorithms for agglomerative hierarchical clustering methods"
133  *
134  * What we found in the past:
135  * - We seem to hit the O(N^3)/worst case of Alg 1
136  * quite often.
137  * - Alg 2 had significant overhead
138  * from the required data-structures.
139  * We could revisit these in the future.
140  */
141 
142 /*
143  * Pairwise distances implementation:
144  * Assuming N total elements 0... N-1,
145  * the pairwise distance matrix
146  * can be represented in a triangular array: <br>
147  * [ (1,0) ] <br>
148  * [ (2,0), (2,1) ] <br>
149  * [ (3,0), (3,1), (3,2)] <br>
150  * [ (4,0), (4,1), (4,2) , (4,3) <br>
151  * [.............................] <br>
152  * [(N-1,0),(N-1,1),(N-1,2),(N-1,3) ... (N-1,N-2)]<br>
153  * With size 1+2+3+ .... (N-1) = N*(N-1)/2
154  *
155  * The lexicographical storage allocation function is
156  * Loc(i,j) = i*(i-1)/2 + j <br>
157  * e.g : <br>
158  * (1,0) => 1 *(1-1)/2 + 0 => 0 <br>
159  * (2,0) => 2 *(2-1)/2 + 0 => 1 <br>
160  * (2,1) => 2 *(2-1)/2 + 1 => 2 <br>
161  * (3,0) => 3 * (3-1)/2 +0 => 3 <br>
162  * Leading to <br>
163  * [(1,0),(2,0),(2,1),(3,0),(3,1),(3,2).... (N-1,N-2)]
164  *
165  * The N-1 Rows map to the value K of the 1st element in the pair
166  * 1,2,3,..,N-1. <br>
167  * Each Row has size K and starts at array positions K*(K-1)/2 <br>
168  * e.g <br>
169  * The row for element 1 starts at array position 0. <br>
170  * The row for element 2 starts at array position 1. <br>
171  * The row for element N-1 starts at array positon (N-1)*(N-2)/2 <br>
172  *
173  * The N-1 Columns map to the value K of the second element in the pair <br>
174  * K= 0,1,2 .., N-2 <br>
175  * The array positions follows (i-1)*i/2+K <br>
176  * where i : K+1 .... N-1 [for(i=K+1;i<N;++i) <br>
177  * e.g <br>
178  * 0 appears as 2nd element in the pair at array positions [0,1,3,6...] <br>
179  * 1 appears as 2nd element in the pair at array positions [2,4,7...] <br>
180  * 2 appears as 2nd element in the pair at array positions [5,8,12....] <br>
181  * N-2 appears as 2nd element once at position [N(N-1)/2-1] <br>
182  */
183 
188 struct triangularToIJ
189 {
190  int I = -1;
191  int J = -1;
192 };
199 inline triangularToIJ
200 convert(int idx)
201 {
202  if (idx<0){
203  throw std::out_of_range("KLGaussianMixtureReduction.cxx::convert : idx is negative");
204  }
205  int i = std::floor((std::sqrt(1 + 8 * idx) + 1) / 2);
206  int j = idx - (i - 1) * i / 2;
207  return {i, j};
208 }
209 
216 inline void
217 calculateAllDistances(const Component1D* componentsIn,
218  float* distancesIn,
219  const int n)
220 {
221  const Component1D* components = std::assume_aligned<GSFConstants::alignment>(componentsIn);
222  float* distances = std::assume_aligned<GSFConstants::alignment>(distancesIn);
223  for (int i = 1; i < n; ++i) {
224  const int indexConst = (i - 1) * i / 2;
225  const Component1D componentI = components[i];
226  for (int j = 0; j < i; ++j) {
227  const Component1D componentJ = components[j];
228  distances[indexConst + j] = symmetricKL(componentI, componentJ);
229  }
230  }
231 }
232 
242 inline int
243 updateDistances(
244  Component1D* ATH_RESTRICT componentsIn,
245  int* ATH_RESTRICT mergingIndexIn,
246  float* ATH_RESTRICT distancesIn,
247  int minFrom,
248  int minTo,
249  int n)
250 {
251  float* distances = std::assume_aligned<GSFConstants::alignment>(distancesIn);
252  Component1D* components = std::assume_aligned<GSFConstants::alignment>(componentsIn);
253  int* mergingIndex = std::assume_aligned<GSFConstants::alignment>(mergingIndexIn);
254  // We swap the last elements with the ones indexed by minFrom.
255  // After this the remaining components we care about
256  // are n-1 which we return
257  const int last = (n - 1);
258  const int indexOffsetJ = (minFrom - 1) * minFrom / 2;
259  const int indexOffsetLast = (last - 1) * last / 2;
260  // we do no need to swap the last with itself
261  if (minFrom != last) {
262  // Rows in distance matrix
263  for (int i = 0; i < minFrom; ++i) {
264  std::swap(distances[indexOffsetJ + i], distances[indexOffsetLast + i]);
265  }
266  // Columns in distance matrix
267  for (int i = minFrom + 1; i < last; ++i) {
268  const int index = (i-1) * i/2 + minFrom;
269  std::swap(distances[index], distances[indexOffsetLast + i]);
270  }
271  // swap the components
272  std::swap(components[minFrom], components[last]);
273  std::swap(mergingIndex[minFrom], mergingIndex[last]);
274  }
275  // In case minTo was indexing the last now it should
276  // be indexing minFrom due to the above swapping
277  if (minTo == last) {
278  minTo = minFrom;
279  }
280  const int indexConst = (minTo - 1) * minTo/2;
281  // This is the component that has been updated
282  const Component1D componentJ = components[minTo];
283  // Rows in distance matrix
284  for (int i = 0; i < minTo; ++i) {
285  const Component1D componentI = components[i];
286  const int index = indexConst + i;
287  distances[index] = symmetricKL(componentI, componentJ);
288  }
289  // Columns in distance matrix
290  for (int i = minTo + 1; i < last; ++i) {
291  const Component1D componentI = components[i];
292  const int index = (i- 1) * i/2 + minTo;
293  distances[index] = symmetricKL(componentI, componentJ);
294  }
295  return last;
296 }
297 
298 } // anonymous namespace with implementation
299 
300 namespace GSFUtils {
305 findMerges(Component1DArray&& componentsIn, const int reducedSize)
306 {
307  const int n = componentsIn.size();
308  int nn = n * (n - 1) / 2;
309  int nnpadded = vAlgs::numPadded<STRIDEForKL>(nn);
310  Component1DArray components(std::move(componentsIn));
311  // Based on the inputSize n allocate enough space for the pairwise distances
313  // initial distance calculation
314  calculateAllDistances(components.buffer(), distances.buffer(), n);
315  // As we merge keep track where things moved
317  std::iota(mergingIndex.begin(), mergingIndex.end(), 0);
318  // Result to be returned
319  MergeArray result{};
320  result.reserve(n);
321  int numberOfComponentsLeft = n;
322  // merge loop
323  while (numberOfComponentsLeft > reducedSize) {
324  // find pair with minimum distance
325  const int minIndex = KLReductionFMV::vIdxOfMin(distances.buffer(), nnpadded);
326  const triangularToIJ conversion = convert(minIndex);
327  int minTo = conversion.I;
328  int minFrom = conversion.J;
329  // This is the convention we had so retained.
330  if (mergingIndex[minTo] < mergingIndex[minFrom]) {
331  std::swap(minTo, minFrom);
332  }
333  // prepare what to return
334  const int miniToreturn = mergingIndex[minTo];
335  const int minjToreturn = mergingIndex[minFrom];
336  result.push_back({ miniToreturn, minjToreturn });
337  // Combine
338  combine(components[minTo], components[minFrom]);
339  // update distances
340  numberOfComponentsLeft = updateDistances(components.buffer(),
341  mergingIndex.buffer(),
342  distances.buffer(),
343  minFrom,
344  minTo,
345  numberOfComponentsLeft);
346 
347  // number of remaining distances padded
348  nn = (numberOfComponentsLeft - 1) * numberOfComponentsLeft / 2;
349  nnpadded = numDistances(nn, distances.buffer());
350  } // end of merge while
351  return result;
352 }
353 
354 } // end namespace GSFUtils
get_generator_info.result
result
Definition: get_generator_info.py:21
KLReductionFMV::vIdxOfMin
int vIdxOfMin(const float *distancesIn, int n)
Definition: KLGaussianMixtureReduction.cxx:37
index
Definition: index.py:1
max
constexpr double max()
Definition: ap_fixedTest.cxx:33
GSFUtils::MergeArray
std::vector< Merge > MergeArray
Definition: KLGaussianMixtureReduction.h:47
ATH_RESTRICT
#define ATH_RESTRICT
Definition: restrict.h:31
GSFUtils::AlignedDynArray::buffer
pointer buffer() noexcept
Get the underlying buffer.
lumiFormat.i
int i
Definition: lumiFormat.py:85
beamspotman.n
n
Definition: beamspotman.py:729
KLReductionFMV
Definition: KLGaussianMixtureReduction.cxx:25
WriteCalibToCool.swap
swap
Definition: WriteCalibToCool.py:94
restrict.h
Macro wrapping the nonstandard restrict keyword.
GSFUtils::AlignedDynArray::begin
iterator begin() noexcept
iterator pointing to the first element
GSFUtils::AlignedDynArray::end
iterator end() noexcept
iterator pointing to the past-the-end element
TMVAToMVAUtils::convert
std::unique_ptr< MVAUtils::BDT > convert(TMVA::MethodBDT *bdt, bool isRegression=true, bool useYesNoLeaf=false)
Definition: TMVAToMVAUtils.h:114
GSFUtils::findMerges
MergeArray findMerges(Component1DArray &&componentsIn, const int reducedSize)
Return which components need to be merged.
Definition: KLGaussianMixtureReduction.cxx:305
KLGaussianMixtureReduction.h
Utilities to facilitate the calculation of the KL divergence/distance between components of the mixtu...
DeMoScan.index
string index
Definition: DeMoScan.py:362
vec.h
Vectorization helpers.
lumiFormat.fill
fill
Definition: lumiFormat.py:104
copySelective.target
string target
Definition: copySelective.py:36
GsfFindIndexOfMinimum.h
LArNewCalib_DelayDump_OFC_Cali.idx
idx
Definition: LArNewCalib_DelayDump_OFC_Cali.py:69
GSFUtils
Dynamic array fullfilling alignment requirements.
Definition: KLGaussianMixtureReduction.cxx:300
I
#define I(x, y, z)
Definition: MD5.cxx:116
GSFUtils::Component1D
struct representing 1D component
Definition: KLGaussianMixtureReduction.h:31
FlavorTagInference::combine
size_t combine(size_t lhs, size_t rhs)
Definition: hash.h:21
GSFUtils::AlignedDynArray
A wrapper around std::aligned_alloc.
Definition: AlignedDynArray.h:30