53 const double meanDifference = componentI.mean - componentJ.mean;
54 const double inverCovSum = componentI.invCov + componentJ.invCov;
55 const double term1 = componentI.invCov * componentJ.cov;
56 const double term2 = componentJ.invCov * componentI.cov;
57 const double term3 = meanDifference * inverCovSum * meanDifference;
58 return term1 + term2 + term3;
72 const double sumWeight = updated.weight + removed.weight;
74 const double invSumWeight = 1. / sumWeight;
75 const double weightI_IJ = updated.weight * invSumWeight;
76 const double weightJ_IJ = removed.weight * invSumWeight;
77 const double meanDiff = (updated.mean - removed.mean);
79 const double sumMean = weightI_IJ * updated.mean + weightJ_IJ * removed.mean;
81 const double sumVariance = weightI_IJ * updated.cov +
82 weightJ_IJ * removed.cov +
83 weightI_IJ * weightJ_IJ * meanDiff * meanDiff;
84 updated.mean = sumMean;
85 updated.cov = sumVariance;
86 updated.invCov = 1. / sumVariance;
87 updated.weight = sumWeight;
95 constexpr
inline int32_t
96 numPadded(
const int32_t
n)
105 numDistances(
const int32_t
n,
float* distancesIn)
107 const int32_t npadded = numPadded(
n);
190 constexpr std::array<int32_t, GSFConstants::maxComponentsAfterConvolution>
193 std::array<int32_t, n>
tmp = {};
194 for (int32_t
i = 0;
i <
n; ++
i) {
204 struct triangularToIJ
214 inline triangularToIJ
223 static const std::vector<triangularToIJ> preMap = []() {
225 constexpr
size_t nn =
n * (
n - 1) / 2;
226 std::vector<triangularToIJ> indexMap(nn);
227 for (int8_t
i = 1;
i <
n; ++
i) {
228 const int32_t indexConst =
offset[
i];
229 for (int8_t j = 0; j <
i; ++j) {
230 indexMap[indexConst + j] = {
i, j };
245 calculateAllDistances(
const Component1D* componentsIn,
250 std::assume_aligned<GSFConstants::alignment>(componentsIn);
252 std::assume_aligned<GSFConstants::alignment>(distancesIn);
253 for (int32_t
i = 1;
i <
n; ++
i) {
254 const int32_t indexConst =
offset[
i];
256 for (int32_t j = 0; j <
i; ++j) {
258 distances[indexConst + j] = symmetricKL(componentI, componentJ);
275 std::array<int8_t, GSFConstants::maxComponentsAfterConvolution>& mergingIndex,
282 std::assume_aligned<GSFConstants::alignment>(distancesIn);
284 std::assume_aligned<GSFConstants::alignment>(componentsIn);
288 const int32_t last = (
n - 1);
289 const int32_t indexOffsetJ =
offset[minFrom];
290 const int32_t indexOffsetLast =
offset[last];
292 if (minFrom != last) {
294 for (int32_t
i = 0;
i < minFrom; ++
i) {
295 std::swap(distances[indexOffsetJ +
i], distances[indexOffsetLast +
i]);
298 for (int32_t
i = minFrom + 1;
i < last; ++
i) {
303 std::swap(components[minFrom], components[last]);
304 std::swap(mergingIndex[minFrom], mergingIndex[last]);
311 const int32_t indexConst =
offset[minTo];
315 for (int32_t
i = 0;
i < minTo; ++
i) {
317 const int32_t
index = indexConst +
i;
318 distances[
index] = symmetricKL(componentI, componentJ);
321 for (int32_t
i = minTo + 1;
i < last; ++
i) {
324 distances[
index] = symmetricKL(componentI, componentJ);
335 const int8_t reducedSize)
339 Component1D* components = std::assume_aligned<GSFConstants::alignment>(
340 copyComponents.components.data());
342 int32_t nn =
n * (
n - 1) / 2;
343 int32_t nnpadded = numPadded(nn);
347 calculateAllDistances(components, distances.buffer(),
n);
349 std::array<int8_t, GSFConstants::maxComponentsAfterConvolution>
351 std::iota(mergingIndex.begin(), mergingIndex.end(), 0);
354 int32_t numberOfComponentsLeft =
n;
356 while (numberOfComponentsLeft > reducedSize) {
358 const int32_t minIndex =
359 findIdxOfMinimum::impl<findIdxOfMinimum::VecMinThenIdx>(
360 distances.buffer(), nnpadded);
361 const triangularToIJ conversion =
convert(minIndex);
362 int8_t minTo = conversion.I;
363 int8_t minFrom = conversion.J;
365 if (mergingIndex[minTo] < mergingIndex[minFrom]) {
369 const int8_t miniToreturn = mergingIndex[minTo];
370 const int8_t minjToreturn = mergingIndex[minFrom];
371 result.merges[
result.numMerges] = { miniToreturn, minjToreturn };
374 combine(components[minTo], components[minFrom]);
376 numberOfComponentsLeft = updateDistances(components,
381 numberOfComponentsLeft);
384 nn =
offset[numberOfComponentsLeft];
385 nnpadded = numDistances(nn, distances.buffer());
399 throw std::runtime_error(
"findMerges :Invalid InputSize or reducedSize");
401 return findMergesImpl(componentsIn,
n, reducedSize);