Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:11:00

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