Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-07-11 07:49:54

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008 
0009 #pragma once
0010 
0011 #include "Acts/Seeding/CandidatesForMiddleSp.hpp"
0012 
0013 #include <limits>
0014 
0015 namespace Acts {
0016 
0017 template <SatisfyCandidateConcept external_space_point_t>
0018 inline void CandidatesForMiddleSp<external_space_point_t>::setMaxElements(
0019     std::size_t nLow, std::size_t nHigh) {
0020   m_maxSizeHigh = nHigh;
0021   m_maxSizeLow = nLow;
0022 
0023   // protection against default numbers
0024   // it may cause std::bad_alloc if we don't protect
0025   if (nHigh == std::numeric_limits<std::size_t>::max() ||
0026       nLow == std::numeric_limits<std::size_t>::max()) {
0027     return;
0028   }
0029 
0030   // Reserve enough memory for all collections
0031   m_storage.reserve(nLow + nHigh);
0032   m_indicesHigh.reserve(nHigh);
0033   m_indicesLow.reserve(nLow);
0034 }
0035 
0036 template <SatisfyCandidateConcept external_space_point_t>
0037 inline void CandidatesForMiddleSp<external_space_point_t>::pop(
0038     std::vector<std::size_t>& indices, std::size_t& currentSize) {
0039   // Remove the candidate with the lowest weight in the collection
0040   // By construction, this element is always the first element in the
0041   // collection.
0042   // The removal works this way: the first element of the collection
0043   // is overwritten with the last element of the collection.
0044   // The number of stored elements (currentSize) is lowered by 1
0045   // The new first element is moved down in the tree to its correct position
0046   std::swap(indices[0], indices[currentSize - 1]);
0047   bubbledw(indices, 0, --currentSize);
0048 }
0049 
0050 template <SatisfyCandidateConcept external_space_point_t>
0051 inline bool CandidatesForMiddleSp<external_space_point_t>::exists(
0052     const std::size_t n, const std::size_t maxSize) const {
0053   // If the element exists, its index is lower than the current number
0054   // of stored elements
0055   return n < maxSize;
0056 }
0057 
0058 template <SatisfyCandidateConcept external_space_point_t>
0059 inline float CandidatesForMiddleSp<external_space_point_t>::weight(
0060     const std::vector<std::size_t>& indices, std::size_t n) const {
0061   // Get the weight of the n-th element
0062   return m_storage[indices[n]].weight;
0063 }
0064 
0065 template <SatisfyCandidateConcept external_space_point_t>
0066 inline void CandidatesForMiddleSp<external_space_point_t>::clear() {
0067   // do not clear max size, this is set only once
0068   // reset to 0 the number of stored elements
0069   m_nHigh = 0;
0070   m_nLow = 0;
0071   // Reset all values to default
0072   m_storage.clear();
0073   m_indicesHigh.clear();
0074   m_indicesLow.clear();
0075 }
0076 
0077 template <SatisfyCandidateConcept external_space_point_t>
0078 bool CandidatesForMiddleSp<external_space_point_t>::push(
0079     external_space_point_t& spB, external_space_point_t& spM,
0080     external_space_point_t& spT, float weight, float zOrigin, bool isQuality) {
0081   // Decide in which collection this candidate may be added to according to the
0082   // isQuality boolean
0083   if (isQuality) {
0084     return push(m_indicesHigh, m_nHigh, m_maxSizeHigh, spB, spM, spT, weight,
0085                 zOrigin, isQuality);
0086   }
0087   return push(m_indicesLow, m_nLow, m_maxSizeLow, spB, spM, spT, weight,
0088               zOrigin, isQuality);
0089 }
0090 
0091 template <SatisfyCandidateConcept external_space_point_t>
0092 bool CandidatesForMiddleSp<external_space_point_t>::push(
0093     std::vector<std::size_t>& indices, std::size_t& n, const std::size_t nMax,
0094     external_space_point_t& spB, external_space_point_t& spM,
0095     external_space_point_t& spT, float weight, float zOrigin, bool isQuality) {
0096   // If we do not want to store candidates, returns
0097   if (nMax == 0) {
0098     return false;
0099   }
0100 
0101   // if there is still space, add anything
0102   if (n < nMax) {
0103     addToCollection(indices, n, nMax,
0104                     value_type(spB, spM, spT, weight, zOrigin, isQuality));
0105     return true;
0106   }
0107 
0108   // if no space, replace one if quality is enough
0109   // compare to element with lowest weight
0110   if (float lowestWeight = this->weight(indices, 0); weight <= lowestWeight) {
0111     return false;
0112   }
0113 
0114   // remove element with lower weight and add this one
0115   pop(indices, n);
0116   addToCollection(indices, n, nMax,
0117                   value_type(spB, spM, spT, weight, zOrigin, isQuality));
0118   return true;
0119 }
0120 
0121 template <SatisfyCandidateConcept external_space_point_t>
0122 void CandidatesForMiddleSp<external_space_point_t>::addToCollection(
0123     std::vector<std::size_t>& indices, std::size_t& n, const std::size_t nMax,
0124     value_type&& element) {
0125   // adds elements to the end of the collection
0126   if (indices.size() == nMax) {
0127     m_storage[indices[n]] = std::move(element);
0128   } else {
0129     m_storage.push_back(std::move(element));
0130     indices.push_back(m_storage.size() - 1);
0131   }
0132   // Move the added element up in the tree to its correct position
0133   bubbleup(indices, n++);
0134 }
0135 
0136 template <SatisfyCandidateConcept external_space_point_t>
0137 void CandidatesForMiddleSp<external_space_point_t>::bubbledw(
0138     std::vector<std::size_t>& indices, std::size_t n, std::size_t actualSize) {
0139   while (n < actualSize) {
0140     // The collection of indexes are sorted as min heap trees
0141     // left child : 2 * n + 1
0142     // right child: 2 * n + 2
0143     float current = weight(indices, n);
0144     std::size_t leftChild = 2 * n + 1;
0145     std::size_t rightChild = 2 * n + 2;
0146 
0147     // We have to move the current node down the tree to its correct position.
0148     // This is done by comparing its weight with the weights of its two
0149     // children. Few things can happen:
0150     //   - there are no children
0151     //   - the current weight is lower than the weight of the children
0152     //   - at least one of the children has a lower weight
0153     // In the first two cases we stop, since we are already in the correct
0154     // position
0155 
0156     // if there is no left child, that also means no right child is present.
0157     // We do nothing
0158     if (!exists(leftChild, actualSize)) {
0159       break;
0160     }
0161 
0162     // At least one of the child is present. Left child for sure, right child we
0163     // have to check. We take the lowest weight of the children. By default this
0164     // is the weight of the left child, and we then check for the right child
0165 
0166     float weightLeftChild = weight(indices, leftChild);
0167 
0168     std::size_t selectedChild = leftChild;
0169     float selectedWeight = weightLeftChild;
0170 
0171     // Check which child has the lower weight
0172     if (exists(rightChild, actualSize)) {
0173       float weightRightChild = weight(indices, rightChild);
0174       if (weightRightChild <= weightLeftChild) {
0175         selectedChild = rightChild;
0176         selectedWeight = weightRightChild;
0177       }
0178     }
0179 
0180     // At this point we have the minimum weight of the children
0181     // We can compare this to the current weight
0182     // If weight of the children is higher we stop
0183     if (selectedWeight >= current) {
0184       break;
0185     }
0186 
0187     // swap and repeat the process
0188     std::swap(indices[n], indices[selectedChild]);
0189     n = selectedChild;
0190   }  // while loop
0191 }
0192 
0193 template <SatisfyCandidateConcept external_space_point_t>
0194 void CandidatesForMiddleSp<external_space_point_t>::bubbleup(
0195     std::vector<std::size_t>& indices, std::size_t n) {
0196   while (n != 0) {
0197     // The collection of indexes are sorted as min heap trees
0198     // parent: (n - 1) / 2;
0199     // this works because it is an integer operation
0200     std::size_t parentIdx = (n - 1) / 2;
0201 
0202     float weightCurrent = weight(indices, n);
0203     // If weight of the parent is lower than this one, we stop
0204     if (float weightParent = weight(indices, parentIdx);
0205         weightParent <= weightCurrent) {
0206       break;
0207     }
0208 
0209     // swap and repeat the process
0210     std::swap(indices[n], indices[parentIdx]);
0211     n = parentIdx;
0212   }
0213 }
0214 
0215 template <SatisfyCandidateConcept external_space_point_t>
0216 std::vector<typename CandidatesForMiddleSp<external_space_point_t>::value_type>
0217 CandidatesForMiddleSp<external_space_point_t>::storage() {
0218   // this will retrieve the entire storage
0219   // the resulting vector is already sorted from high to low quality
0220   std::vector<value_type> output(m_nHigh + m_nLow);
0221   std::size_t outIdx = output.size() - 1;
0222 
0223   // rely on the fact that m_indices* are both min heap trees
0224   // Sorting comes naturally by popping elements one by one and
0225   // placing this element at the end of the output vector
0226   while (m_nHigh != 0 || m_nLow != 0) {
0227     // no entries in collection high, we attach the entire low collection
0228     if (m_nHigh == 0) {
0229       std::size_t idx = m_nLow;
0230       for (std::size_t i(0); i < idx; i++) {
0231         output[outIdx--] = std::move(m_storage[m_indicesLow[0]]);
0232         pop(m_indicesLow, m_nLow);
0233       }
0234       break;
0235     }
0236 
0237     // no entries in collection low, we attach the entire high collection
0238     if (m_nLow == 0) {
0239       std::size_t idx = m_nHigh;
0240       for (std::size_t i(0); i < idx; i++) {
0241         output[outIdx--] = std::move(m_storage[m_indicesHigh[0]]);
0242         pop(m_indicesHigh, m_nHigh);
0243       }
0244       break;
0245     }
0246 
0247     // Both have entries, get the minimum
0248     if (descendingByQuality(m_storage[m_indicesLow[0]],
0249                             m_storage[m_indicesHigh[0]])) {
0250       output[outIdx--] = std::move(m_storage[m_indicesHigh[0]]);
0251       pop(m_indicesHigh, m_nHigh);
0252     } else {
0253       output[outIdx--] = std::move(m_storage[m_indicesLow[0]]);
0254       pop(m_indicesLow, m_nLow);
0255     }
0256 
0257   }  // while loop
0258 
0259   clear();
0260   return output;
0261 }
0262 
0263 template <SatisfyCandidateConcept external_space_point_t>
0264 bool CandidatesForMiddleSp<external_space_point_t>::descendingByQuality(
0265     const value_type& i1, const value_type& i2) {
0266   if (i1.weight != i2.weight) {
0267     return i1.weight > i2.weight;
0268   }
0269 
0270   // This is for the case when the weights from different seeds
0271   // are same. This makes cpu & cuda results same
0272 
0273   const auto& bottomL1 = i1.bottom;
0274   const auto& middleL1 = i1.middle;
0275   const auto& topL1 = i1.top;
0276 
0277   const auto& bottomL2 = i2.bottom;
0278   const auto& middleL2 = i2.middle;
0279   const auto& topL2 = i2.top;
0280 
0281   float seed1_sum = 0.;
0282   float seed2_sum = 0.;
0283 
0284   seed1_sum += bottomL1->y() * bottomL1->y() + bottomL1->z() * bottomL1->z();
0285   seed1_sum += middleL1->y() * middleL1->y() + middleL1->z() * middleL1->z();
0286   seed1_sum += topL1->y() * topL1->y() + topL1->z() * topL1->z();
0287 
0288   seed2_sum += bottomL2->y() * bottomL2->y() + bottomL2->z() * bottomL2->z();
0289   seed2_sum += middleL2->y() * middleL2->y() + middleL2->z() * middleL2->z();
0290   seed2_sum += topL2->y() * topL2->y() + topL2->z() * topL2->z();
0291 
0292   return seed1_sum > seed2_sum;
0293 }
0294 
0295 template <SatisfyCandidateConcept external_space_point_t>
0296 bool CandidatesForMiddleSp<external_space_point_t>::ascendingByQuality(
0297     const value_type& i1, const value_type& i2) {
0298   return !descendingByQuality(i1, i2);
0299 }
0300 
0301 template <SatisfyCandidateConcept external_space_point_t>
0302 std::size_t
0303 CandidatesForMiddleSp<external_space_point_t>::nLowQualityCandidates() const {
0304   return m_nLow;
0305 }
0306 
0307 template <SatisfyCandidateConcept external_space_point_t>
0308 std::size_t
0309 CandidatesForMiddleSp<external_space_point_t>::nHighQualityCandidates() const {
0310   return m_nHigh;
0311 }
0312 
0313 }  // namespace Acts