VecGeom/navigation/HybridSafetyEstimator.h

0001 /*
0002  *  HybridSafetyEstimator.h
0003  *
0004  *  Created on: 22.11.2015
0005  *      Author: sandro.wenzel@cern.ch
0006  *
0007  *  (based on prototype implementation by Yang Zhang (Sep 2015)
0008  */
0009
0010 #ifndef NAVIGATION_HYBRIDSAFETYESTIMATOR_H_
0011 #define NAVIGATION_HYBRIDSAFETYESTIMATOR_H_
0012
0013 #include "VecGeom/navigation/VSafetyEstimator.h"
0014 #include "VecGeom/management/HybridManager2.h"
0015
0016 namespace vecgeom {
0017 inline namespace VECGEOM_IMPL_NAMESPACE {
0018
0019 //! a safety estimator using a (vectorized) search through bounding boxes to exclude certain daughter volumes
0020 //! to talk to
0021 class HybridSafetyEstimator : public VSafetyEstimatorHelper<HybridSafetyEstimator> {
0022
0023 private:
0024   // we keep a reference to the ABBoxManager ( avoids calling Instance() on this guy all the time )
0025   HybridManager2 &fAccelerationStructureManager;
0026
0027   HybridSafetyEstimator()
0028       : VSafetyEstimatorHelper<HybridSafetyEstimator>(), fAccelerationStructureManager(HybridManager2::Instance())
0029   {
0030   }
0031
0032   // convert index to physical daugher
0033   VPlacedVolume const *LookupDaughter(LogicalVolume const *lvol, int id) const
0034   {
0035     VECGEOM_VALIDATE(id >= 0, << "access with negative index");
0036     VECGEOM_VALIDATE(size_t(id) < lvol->GetDaughtersp()->size(), << "access beyond size of daughterlist ");
0037     return lvol->GetDaughtersp()->operator[](id);
0038   }
0039
0040 public:
0041   // a simple sort class (based on insertionsort)
0042   template <typename T> //, typename Cmp>
0043   static void insertionsort(T *arr, unsigned int N)
0044   {
0045     for (unsigned short i = 1; i < N; ++i) {
0046       T value    = arr[i];
0047       short hole = i;
0048
0049       for (; hole > 0 && value.second < arr[hole - 1].second; --hole)
0050         arr[hole] = arr[hole - 1];
0051
0052       arr[hole] = value;
0053     }
0054   }
0055
0056   // helper structure to find the candidate set for safety calculations
0057   size_t GetSafetyCandidates_v(HybridManager2::HybridBoxAccelerationStructure const &accstructure,
0058                                Vector3D<Precision> const &point, HybridManager2::BoxIdDistancePair_t *boxsafetypairs,
0059                                Precision upper_squared_limit) const
0060   {
0061     size_t count = 0;
0062     Vector3D<float> pointfloat((float)point.x(), (float)point.y(), (float)point.z());
0063     int halfvectorsize, numberOfNodes;
0064     auto boxes_v = fAccelerationStructureManager.GetABBoxes_v(accstructure, halfvectorsize, numberOfNodes);
0065     auto const *nodeToDaughters = accstructure.fNodeToDaughters;
0066     constexpr auto kVS          = vecCore::VectorSize<HybridManager2::Float_v>();
0067
0068     for (int index = 0, nodeindex = 0; index < halfvectorsize * 2; index += 2 * (kVS + 1), nodeindex += kVS) {
0069       HybridManager2::Float_v safetytoboxsqr =
0070           ABBoxImplementation::ABBoxSafetySqr(boxes_v[index], boxes_v[index + 1], pointfloat);
0071       auto closer = safetytoboxsqr < HybridManager2::Float_v(upper_squared_limit);
0072       if (!vecCore::MaskEmpty(closer)) {
0073         for (size_t i = 0 /*closer.firstOne()*/; i < kVS; ++i) {
0074           if (vecCore::MaskLaneAt(closer, i)) {
0075             safetytoboxsqr =
0076                 ABBoxImplementation::ABBoxSafetySqr(boxes_v[index + 2 * i + 2], boxes_v[index + 2 * i + 3], pointfloat);
0077             auto closer1 = safetytoboxsqr < HybridManager2::Float_v(upper_squared_limit);
0078             if (!vecCore::MaskEmpty(closer1)) {
0079               for (size_t j = 0 /*closer.firstOne()*/; j < kVS; ++j) { // leaf node
0080                 if (vecCore::MaskLaneAt(closer1, j)) {
0081                   VECGEOM_ASSERT(count < VECGEOM_MAXFACETS);
0082                   boxsafetypairs[count] = HybridManager2::BoxIdDistancePair_t(nodeToDaughters[nodeindex + i][j],
0083                                                                               vecCore::LaneAt(safetytoboxsqr, j));
0084                   count++;
0085                 }
0086               }
0087             }
0088           }
0089         }
0090       }
0091     }
0092     return count;
0093   }
0094
0095   // Improved safety estimator that can dynamically reduce the upper search limit.
0096   // Internally calls ABBoxSafetyRangeSqr that besides validating a candidate
0097   // if closer than the limit, it also gives the upper limit for the range where
0098   // the candidate can be surely found. This allows to update the search limit after
0099   // each check, giving much less final candidates.
0100   //
0101   //         current checked
0102   //         box  _________
0103   // point       |         |          upper_limit
0104   //   x---------|candidate|----------|
0105   //             |_________|
0106   //   |-------------------|
0107   //                       updated upper_limit
0108
0109   size_t GetSafetyCandidates2_v(HybridManager2::HybridBoxAccelerationStructure const &accstructure,
0110                                 Vector3D<Precision> const &point, HybridManager2::BoxIdDistancePair_t *hitlist,
0111                                 Precision upper_squared_limit) const
0112   {
0113     using Float_v = vecgeom::VectorBackend::Float_v;
0114     size_t count  = 0;
0115     int numberOfNodes, size;
0116     auto boxes_v                = fAccelerationStructureManager.GetABBoxes_v(accstructure, size, numberOfNodes);
0117     constexpr auto kVS          = vecCore::VectorSize<HybridManager2::Float_v>();
0118     auto const *nodeToDaughters = accstructure.fNodeToDaughters;
0119
0120     Vector3D<float> pointfloat((float)point.x(), (float)point.y(), (float)point.z());
0121
0122     // Loop internal nodes (internal node = kVS clusters)
0123     for (size_t index = 0, nodeindex = 0; index < size_t(size) * 2; index += 2 * (kVS + 1), nodeindex += kVS) {
0124       // index = start index for internal node
0125       // nodeindex = cluster index
0126       Float_v distmaxsqr; // Maximum distance that may still touch the box
0127       Float_v safetytonodesqr =
0128           ABBoxImplementation::ABBoxSafetyRangeSqr(boxes_v[index], boxes_v[index + 1], pointfloat, distmaxsqr);
0129       // Find minimum of distmaxsqr
0130       for (size_t i = 0; i < kVS; ++i) {
0131         Precision distmaxsqr_s = vecCore::LaneAt(distmaxsqr, i);
0132         if (distmaxsqr_s < upper_squared_limit) upper_squared_limit = distmaxsqr_s;
0133       }
0134       auto hit = safetytonodesqr < ABBoxManager<Precision>::Real_s(upper_squared_limit);
0135       if (!vecCore::MaskEmpty(hit)) {
0136         for (size_t i = 0; i < kVS; ++i) {
0137           if (vecCore::MaskLaneAt(hit, i)) {
0138             Float_v safetytoboxsqr = ABBoxImplementation::ABBoxSafetyRangeSqr(
0139                 boxes_v[index + 2 * (i + 1)], boxes_v[index + 2 * (i + 1) + 1], pointfloat, distmaxsqr);
0140
0141             auto hit1 = safetytoboxsqr < ABBoxManager<Precision>::Real_s(upper_squared_limit);
0142             if (!vecCore::MaskEmpty(hit1)) {
0143               // loop bounding boxes in the cluster
0144               for (size_t j = 0; j < kVS; ++j) {
0145                 if (vecCore::MaskLaneAt(hit1, j)) {
0146                   VECGEOM_ASSERT(count < VECGEOM_MAXFACETS);
0147                   hitlist[count]         = HybridManager2::BoxIdDistancePair_t(nodeToDaughters[nodeindex + i][j],
0148                                                                                vecCore::LaneAt(safetytoboxsqr, j));
0149                   Precision distmaxsqr_s = vecCore::LaneAt(distmaxsqr, j);
0150                   // Reduce the upper limit
0151                   if (distmaxsqr_s < upper_squared_limit) upper_squared_limit = distmaxsqr_s;
0152                   count++;
0153                 }
0154               }
0155             }
0156           }
0157         }
0158       }
0159     }
0160     return count;
0161   }
0162
0163   template <typename AccStructure, typename Func>
0164   VECGEOM_FORCE_INLINE void BVHSortedSafetyLooper(AccStructure const &accstructure,
0165                                                   Vector3D<Precision> const &localpoint, Func &&userhook,
0166                                                   Precision upper_squared_limit) const
0167   {
0168     // The following construct reserves stackspace for objects
0169     // of type IdDistPair_t WITHOUT initializing those objects
0170     using IdDistPair_t = HybridManager2::BoxIdDistancePair_t;
0171     char stackspace[VECGEOM_MAXFACETS * sizeof(IdDistPair_t)];
0172     IdDistPair_t *hitlist = reinterpret_cast<IdDistPair_t *>(&stackspace);
0173
0174     // Get candidates using HybridSafetyEstimator
0175     //    HybridSafetyEstimator *hse = static_cast<HybridSafetyEstimator*>(HybridSafetyEstimator::Instance());
0176     //    auto ncandidates = hse->GetSafetyCandidates_v(accstructure, localpoint, hitlist, upper_squared_limit);
0177
0178     // Get candidates in vectorized mode
0179     auto ncandidates = GetSafetyCandidates2_v(accstructure, localpoint, hitlist, upper_squared_limit);
0180
0181     // sort candidates according to their bounding volume safety distance
0182     insertionsort(hitlist, ncandidates);
0183
0184     for (size_t index = 0; index < ncandidates; ++index) {
0185       auto hitbox = hitlist[index];
0186       // here we got the hit candidates
0187       // now we execute user specific code to process this "hitbox"
0188       auto done = userhook(hitbox);
0189       if (done) break;
0190     }
0191   }
0192
0193   static constexpr const char *gClassNameString = "HybridSafetyEstimator";
0194
0195   VECGEOM_FORCE_INLINE
0196   VECCORE_ATT_HOST_DEVICE
0197   Precision TreatSafetyToIn(Vector3D<Precision> const &localpoint, LogicalVolume const *lvol, Precision outsafety) const
0198   {
0199     // a stack based workspace array
0200     // The following construct reserves stackspace for objects
0201     // of type IdDistPair_t WITHOUT initializing those objects
0202     using IdDistPair_t = HybridManager2::BoxIdDistancePair_t;
0203     char unitstackspace[VECGEOM_MAXDAUGHTERS * sizeof(IdDistPair_t)];
0204     IdDistPair_t *boxsafetylist = reinterpret_cast<IdDistPair_t *>(&unitstackspace);
0205
0206     double safety    = outsafety; // we use the outsafety estimate as starting point
0207     double safetysqr = safety * safety;
0208
0209     // safety to bounding boxes
0210     if (safety > 0. && lvol->GetDaughtersp()->size() > 0) {
0211       // calculate squared bounding box safeties in vectorized way
0212       auto ncandidates = GetSafetyCandidates_v(*fAccelerationStructureManager.GetAccStructure(lvol), localpoint,
0213                                                boxsafetylist, safetysqr);
0214       // not sorting the candidate list ( which one could do )
0215       for (unsigned int candidate = 0; candidate < ncandidates; ++candidate) {
0216         auto boxsafetypair = boxsafetylist[candidate];
0217         if (boxsafetypair.second < safetysqr) {
0218           VPlacedVolume const *cand = LookupDaughter(lvol, boxsafetypair.first);
0219           if (size_t(boxsafetypair.first) > lvol->GetDaughtersp()->size()) break;
0220           auto candidatesafety = cand->SafetyToIn(localpoint);
0221 #ifdef VERBOSE
0222           if (candidatesafety * candidatesafety > boxsafetypair.second && boxsafetypair.second > 0)
0223             std::cerr << "real safety smaller than boxsafety \n";
0224 #endif
0225           if (candidatesafety < safety) {
0226             safety    = candidatesafety;
0227             safetysqr = safety * safety;
0228           }
0229         }
0230       }
0231     }
0232     return safety;
0233   }
0234
0235   // this is (almost) the same code as in SimpleABBoxSafetyEstimator --> avoid this
0236   VECGEOM_FORCE_INLINE
0237   VECCORE_ATT_HOST_DEVICE
0238   virtual Precision ComputeSafetyForLocalPoint(Vector3D<Precision> const &localpoint,
0239                                                VPlacedVolume const *pvol) const override
0240   {
0241     // safety to mother
0242     double safety = pvol->SafetyToOut(localpoint);
0243     if (safety <= 0.) {
0244       return 0.;
0245     }
0246     return TreatSafetyToIn(localpoint, pvol->GetLogicalVolume(), safety);
0247   }
0248
0249   VECGEOM_FORCE_INLINE
0250   VECCORE_ATT_HOST_DEVICE
0251   virtual Precision ComputeSafetyToDaughtersForLocalPoint(Vector3D<Precision> const &localpoint,
0252                                                           LogicalVolume const *lvol) const override
0253   {
0254     return TreatSafetyToIn(localpoint, lvol, kInfLength);
0255   }
0256
0257   static VSafetyEstimator *Instance()
0258   {
0259     static HybridSafetyEstimator instance;
0260     return &instance;
0261   }
0262
0263 }; // end class
0264 } // namespace VECGEOM_IMPL_NAMESPACE
0265 } // namespace vecgeom
0266
0267 #endif /* NAVIGATION_SIMPLEABBOXSAFETYESTIMATOR_H_ */