ATLAS Offline Software
KLGaussianMixtureReduction.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 
8 //
10 //
11 #include "CxxUtils/restrict.h"
12 #include "CxxUtils/vec.h"
13 //
14 #include <cmath>
15 #include <memory>
16 #include <limits>
17 #include <numeric>
18 #include <stdexcept>
19 #include <vector>
20 
28 namespace KLReductionFMV {
29 //clang FMV needs a namespace :/
30 #if HAVE_FUNCTION_MULTIVERSIONING
31 [[gnu::target("default")]]
32 #endif
33 int vIdxOfMin(const float* distancesIn, int n) {
34  return vAlgs::vIdxOfMin<128>(distancesIn, n);
35 }
36 #if HAVE_FUNCTION_MULTIVERSIONING
37 [[gnu::target("avx2")]]
38 int vIdxOfMin(const float* distancesIn, int n) {
39  return vAlgs::vIdxOfMin<256>(distancesIn, n);
40 }
41 #endif
42 } // namespace KLReductionFMV
43 
44 namespace {
45 //internal implementation methods
46 
47 //We want to be using up to a 256 ISA, these cover also a narrower one
48 constexpr size_t STRIDEForKL = vAlgs::strideOfNumSIMDVec<256,float>(4);
49 constexpr size_t ALIGNMENTForKL = vAlgs::alignmentForArray<256>();
50 using namespace GSFUtils;
51 
64 inline float
65 symmetricKL(const Component1D& ATH_RESTRICT componentI,
66  const Component1D& ATH_RESTRICT componentJ)
67 {
68  const double meanDifference = componentI.mean - componentJ.mean;
69  const double inverCovSum = componentI.invCov + componentJ.invCov;
70  const double term1 = componentI.invCov * componentJ.cov;
71  const double term2 = componentJ.invCov * componentI.cov;
72  const double term3 = meanDifference * inverCovSum * meanDifference;
73  return term1 + term2 + term3;
74 }
82 inline void
85 {
86 
87  const double sumWeight = updated.weight + removed.weight;
88 
89  const double invSumWeight = 1. / sumWeight;
90  const double weightI_IJ = updated.weight * invSumWeight;
91  const double weightJ_IJ = removed.weight * invSumWeight;
92  const double meanDiff = (updated.mean - removed.mean);
93 
94  const double sumMean = weightI_IJ * updated.mean + weightJ_IJ * removed.mean;
95 
96  const double sumVariance = weightI_IJ * updated.cov +
97  weightJ_IJ * removed.cov +
98  weightI_IJ * weightJ_IJ * meanDiff * meanDiff;
99  updated.mean = sumMean;
100  updated.cov = sumVariance;
101  updated.invCov = 1. / sumVariance;
102  updated.weight = sumWeight;
103 }
104 
105 inline int
106 numDistances(const int n, float* distancesIn)
107 {
108  const int npadded = vAlgs::numPadded<STRIDEForKL>(n);
109  // Make sure the padded elements are set to max
110  std::fill(
111  distancesIn + n, distancesIn + npadded, std::numeric_limits<float>::max());
112  return npadded;
113 }
114 
115 /*
116  * This is a O(N^3) time N^2 space algorithm.
117  * It can be seen as tailored implementation for
118  * our problem of the basic algorithm
119  * for Hierarchical Clustering.
120  *
121  * We rely on the fast findMinimunIndex above
122  * and a triangular array representation to
123  * reduce the actual time for our problem
124  * (max N=72)
125  *
126  * We opt for fixed size arrays of max N(=72) elements,
127  * but dynamic arrays of ~ N*(N-1)/2 elements.
128  *
129  * Existing alternatives in the literature:
130  * - 1. O(N^3) worst-case time , O(n^2) best
131  * - 2. O(N^2 log(N)) worst-time
132  * See :
133  * "Modern hierarchical, agglomerative clustering algorithms"
134  * https://arxiv.org/abs/1109.2378
135  * or
136  * "Efficient algorithms for agglomerative hierarchical clustering methods"
137  *
138  * What we found in the past:
139  * - We seem to hit the O(N^3)/worst case of Alg 1
140  * quite often.
141  * - Alg 2 had significant overhead
142  * from the required data-structures.
143  * We could revisit these in the future.
144  */
145 
146 /*
147  * Pairwise distances implementation:
148  * Assuming N total elements 0... N-1,
149  * the pairwise distance matrix
150  * can be represented in a triangular array: <br>
151  * [ (1,0) ] <br>
152  * [ (2,0), (2,1) ] <br>
153  * [ (3,0), (3,1), (3,2)] <br>
154  * [ (4,0), (4,1), (4,2) , (4,3) <br>
155  * [.............................] <br>
156  * [(N-1,0),(N-1,1),(N-1,2),(N-1,3) ... (N-1,N-2)]<br>
157  * With size 1+2+3+ .... (N-1) = N*(N-1)/2
158  *
159  * The lexicographical storage allocation function is
160  * Loc(i,j) = i*(i-1)/2 + j <br>
161  * e.g : <br>
162  * (1,0) => 1 *(1-1)/2 + 0 => 0 <br>
163  * (2,0) => 2 *(2-1)/2 + 0 => 1 <br>
164  * (2,1) => 2 *(2-1)/2 + 1 => 2 <br>
165  * (3,0) => 3 * (3-1)/2 +0 => 3 <br>
166  * Leading to <br>
167  * [(1,0),(2,0),(2,1),(3,0),(3,1),(3,2).... (N-1,N-2)]
168  *
169  * The N-1 Rows map to the value K of the 1st element in the pair
170  * 1,2,3,..,N-1. <br>
171  * Each Row has size K and starts at array positions K*(K-1)/2 <br>
172  * e.g <br>
173  * The row for element 1 starts at array position 0. <br>
174  * The row for element 2 starts at array position 1. <br>
175  * The row for element N-1 starts at array positon (N-1)*(N-2)/2 <br>
176  *
177  * The N-1 Columns map to the value K of the second element in the pair <br>
178  * K= 0,1,2 .., N-2 <br>
179  * The array positions follows (i-1)*i/2+K <br>
180  * where i : K+1 .... N-1 [for(i=K+1;i<N;++i) <br>
181  * e.g <br>
182  * 0 appears as 2nd element in the pair at array positions [0,1,3,6...] <br>
183  * 1 appears as 2nd element in the pair at array positions [2,4,7...] <br>
184  * 2 appears as 2nd element in the pair at array positions [5,8,12....] <br>
185  * N-2 appears as 2nd element once at position [N(N-1)/2-1] <br>
186  */
187 
191 constexpr std::array<int, GSFConstants::maxComponentsAfterConvolution>
192  offset = []() {
194  std::array<int, n> tmp = {};
195  for (int i = 0; i < n; ++i) {
196  tmp[i] = (i - 1) * i / 2;
197  }
198  return tmp;
199  }();
200 
205 struct triangularToIJ
206 {
207  int8_t I = -1;
208  int8_t J = -1;
209 };
215 inline triangularToIJ
216 convert(int idx)
217 {
218  // We prefer to preMap the maximum 2556 elements.
219  // Alternatively one can use the following
220  // if pre-mapping becomes an issue
221  // (see https://hal.archives-ouvertes.fr/hal-02047514/document)
222  // int8_t i = std::floor((std::sqrt(1 + 8 * idx) + 1) / 2);
223  // int8_t j = idx - (i - 1) * i / 2;
224  static const std::vector<triangularToIJ> preMap = []() {
226  constexpr size_t nn = n * (n - 1) / 2;
227  std::vector<triangularToIJ> indexMap(nn);
228  for (int8_t i = 1; i < n; ++i) {
229  const int indexConst = offset[i];
230  for (int8_t j = 0; j < i; ++j) {
231  indexMap[indexConst + j] = { i, j };
232  }
233  }
234  return indexMap;
235  }();
236  return preMap[idx];
237 }
238 
245 inline void
246 calculateAllDistances(const Component1D* componentsIn,
247  float* distancesIn,
248  const int n)
249 {
250  const Component1D* components =
251  std::assume_aligned<GSFConstants::alignment>(componentsIn);
252  float* distances =
253  std::assume_aligned<GSFConstants::alignment>(distancesIn);
254  for (int i = 1; i < n; ++i) {
255  const int indexConst = offset[i];
256  const Component1D componentI = components[i];
257  for (int j = 0; j < i; ++j) {
258  const Component1D componentJ = components[j];
259  distances[indexConst + j] = symmetricKL(componentI, componentJ);
260  }
261  }
262 }
263 
273 inline int
274 updateDistances(
275  Component1D* ATH_RESTRICT componentsIn,
276  std::array<int8_t, GSFConstants::maxComponentsAfterConvolution>& mergingIndex,
277  float* ATH_RESTRICT distancesIn,
278  int minFrom,
279  int minTo,
280  int n)
281 {
282  float* distances =
283  std::assume_aligned<GSFConstants::alignment>(distancesIn);
284  Component1D* components =
285  std::assume_aligned<GSFConstants::alignment>(componentsIn);
286  // We swap the last elements with the ones indexed by minFrom.
287  // After this the remaining components we care about
288  // are n-1 which we return
289  const int last = (n - 1);
290  const int indexOffsetJ = offset[minFrom];
291  const int indexOffsetLast = offset[last];
292  // we do no need to swap the last with itself
293  if (minFrom != last) {
294  // Rows in distance matrix
295  for (int i = 0; i < minFrom; ++i) {
296  std::swap(distances[indexOffsetJ + i], distances[indexOffsetLast + i]);
297  }
298  // Columns in distance matrix
299  for (int i = minFrom + 1; i < last; ++i) {
300  const int index = offset[i] + minFrom;
301  std::swap(distances[index], distances[indexOffsetLast + i]);
302  }
303  // swap the components
304  std::swap(components[minFrom], components[last]);
305  std::swap(mergingIndex[minFrom], mergingIndex[last]);
306  }
307  // In case minTo was indexing the last now it should
308  // be indexing minFrom due to the above swapping
309  if (minTo == last) {
310  minTo = minFrom;
311  }
312  const int indexConst = offset[minTo];
313  // This is the component that has been updated
314  const Component1D componentJ = components[minTo];
315  // Rows in distance matrix
316  for (int i = 0; i < minTo; ++i) {
317  const Component1D componentI = components[i];
318  const int index = indexConst + i;
319  distances[index] = symmetricKL(componentI, componentJ);
320  }
321  // Columns in distance matrix
322  for (int i = minTo + 1; i < last; ++i) {
323  const Component1D componentI = components[i];
324  const int index = offset[i] + minTo;
325  distances[index] = symmetricKL(componentI, componentJ);
326  }
327  return last;
328 }
329 
334 findMergesImpl(const Component1DArray& componentsIn,
335  const int n,
336  const int8_t reducedSize)
337 {
338  // copy the array for internal use
339  Component1DArray copyComponents(componentsIn);
340  Component1D* components = std::assume_aligned<GSFConstants::alignment>(
341  copyComponents.components.data());
342  // Based on the inputSize n allocate enough space for the pairwise distances
343  int nn = n * (n - 1) / 2;
344  int nnpadded = vAlgs::numPadded<STRIDEForKL>(nn);
346  nnpadded, std::numeric_limits<float>::max());
347  // initial distance calculation
348  calculateAllDistances(components, distances.buffer(), n);
349  // As we merge keep track where things moved
350  std::array<int8_t, GSFConstants::maxComponentsAfterConvolution>
351  mergingIndex{};
352  std::iota(mergingIndex.begin(), mergingIndex.end(), 0);
353  // Result to be returned
354  MergeArray result{};
355  int numberOfComponentsLeft = n;
356  // merge loop
357  while (numberOfComponentsLeft > reducedSize) {
358  // find pair with minimum distance
359  const int minIndex = KLReductionFMV::vIdxOfMin(distances.buffer(), nnpadded);
360  const triangularToIJ conversion = convert(minIndex);
361  int8_t minTo = conversion.I;
362  int8_t minFrom = conversion.J;
363  // This is the convention we had so retained.
364  if (mergingIndex[minTo] < mergingIndex[minFrom]) {
365  std::swap(minTo, minFrom);
366  }
367  // prepare what to return
368  const int8_t miniToreturn = mergingIndex[minTo];
369  const int8_t minjToreturn = mergingIndex[minFrom];
370  result.merges[result.numMerges] = { miniToreturn, minjToreturn };
371  ++result.numMerges;
372  // Combine
373  combine(components[minTo], components[minFrom]);
374  // update distances
375  numberOfComponentsLeft = updateDistances(components,
376  mergingIndex,
377  distances.buffer(),
378  minFrom,
379  minTo,
380  numberOfComponentsLeft);
381 
382  // number of remaining distances padded
383  nn = offset[numberOfComponentsLeft];
384  nnpadded = numDistances(nn, distances.buffer());
385  } // end of merge while
386  return result;
387 }
388 
389 } // anonymous namespace with implementation
390 
391 namespace GSFUtils {
393 findMerges(const Component1DArray& componentsIn, const int8_t reducedSize)
394 {
395  const int n = componentsIn.numComponents;
397  reducedSize > n) {
398  throw std::runtime_error("findMerges :Invalid InputSize or reducedSize");
399  }
400  return findMergesImpl(componentsIn, n, reducedSize);
401 }
402 } // end namespace GSFUtils
get_generator_info.result
result
Definition: get_generator_info.py:21
GSFUtils::MergeArray
struct representing an array or the merges.
Definition: KLGaussianMixtureReduction.h:54
KLReductionFMV::vIdxOfMin
int vIdxOfMin(const float *distancesIn, int n)
Definition: KLGaussianMixtureReduction.cxx:33
GSFUtils::Component1DArray::numComponents
int numComponents
Definition: KLGaussianMixtureReduction.h:46
index
Definition: index.py:1
max
constexpr double max()
Definition: ap_fixedTest.cxx:33
ATH_RESTRICT
#define ATH_RESTRICT
Definition: restrict.h:31
FlavorTagDiscriminants::combine
size_t combine(size_t lhs, size_t rhs)
Definition: hash.h:21
lumiFormat.i
int i
Definition: lumiFormat.py:85
beamspotman.n
n
Definition: beamspotman.py:731
GSFConstants::maxComponentsAfterConvolution
constexpr int8_t maxComponentsAfterConvolution
The maximum size State x Bethe-Heitler components The typical number we use is the max 6x12 = 72 i....
Definition: GsfConstants.h:61
KLReductionFMV
Definition: KLGaussianMixtureReduction.cxx:28
DeMoUpdate.tmp
string tmp
Definition: DeMoUpdate.py:1167
AlignedDynArray.h
WriteCalibToCool.swap
swap
Definition: WriteCalibToCool.py:94
fill
void fill(H5::Group &out_file, size_t iterations)
Definition: test-hdf5-writer.cxx:95
restrict.h
Macro wrapping the nonstandard restrict keyword.
TMVAToMVAUtils::convert
std::unique_ptr< MVAUtils::BDT > convert(TMVA::MethodBDT *bdt, bool isRegression=true, bool useYesNoLeaf=false)
Definition: TMVAToMVAUtils.h:114
KLGaussianMixtureReduction.h
Utilities to facilitate the calculation of the KL divergence/distance between components of the mixtu...
DeMoScan.index
string index
Definition: DeMoScan.py:364
vec.h
Vectorization helpers.
copySelective.target
string target
Definition: copySelective.py:37
GSFFindIndexOfMinimum.h
convertTimingResiduals.offset
offset
Definition: convertTimingResiduals.py:71
LArNewCalib_DelayDump_OFC_Cali.idx
idx
Definition: LArNewCalib_DelayDump_OFC_Cali.py:69
GSFUtils
Dynamic array fullfilling alignment requirements.
Definition: KLGaussianMixtureReduction.cxx:391
GsfConstants.h
I
#define I(x, y, z)
Definition: MD5.cxx:116
GSFUtils::Component1D
struct representing 1D component
Definition: KLGaussianMixtureReduction.h:30
GSFUtils::findMerges
MergeArray findMerges(const Component1DArray &componentsIn, const int8_t reducedSize)
Find the order in which the components need to be merged.
Definition: KLGaussianMixtureReduction.cxx:393
GSFUtils::Component1DArray
struct representing an array of 1D component.
Definition: KLGaussianMixtureReduction.h:42
GSFUtils::AlignedDynArray
A wrapper around std::aligned_alloc.
Definition: AlignedDynArray.h:30