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     assert(id >= 0 && "access with negative index");
0036     assert(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                   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::Real_t(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::Real_t(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                   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
0165   void BVHSortedSafetyLooper(AccStructure const &accstructure, 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   virtual Real_v ComputeSafetyForLocalPoint(Vector3D<Real_v> const &localpoint, VPlacedVolume const *pvol,
0198                                             Bool_v m) const override
0199   {
0200     using vecCore::AssignLane;
0201     using vecCore::LaneAt;
0202     Real_v safety(0.);
0203     if (!vecCore::MaskEmpty(m)) {
0204       // SIMD safety to mother
0205       safety = pvol->SafetyToOut(localpoint);
0206
0207       LogicalVolume const *lvol = pvol->GetLogicalVolume();
0208       // now loop over the voxelized treatment of safety to in
0209       for (unsigned int i = 0; i < vecCore::VectorSize<Real_v>(); ++i) {
0210         if (vecCore::MaskLaneAt(m, i)) {
0211           AssignLane(safety, i,
0212                      TreatSafetyToIn(Vector3D<Precision>(LaneAt(localpoint.x(), i), LaneAt(localpoint.y(), i),
0213                                                          LaneAt(localpoint.z(), i)),
0214                                      lvol, LaneAt(safety, i)));
0215         } else {
0216           AssignLane(safety, i, 0.);
0217         }
0218       }
0219     }
0220     return safety;
0221   }
0222
0223   VECGEOM_FORCE_INLINE
0224   VECCORE_ATT_HOST_DEVICE
0225   Precision TreatSafetyToIn(Vector3D<Precision> const &localpoint, LogicalVolume const *lvol, Precision outsafety) const
0226   {
0227     // a stack based workspace array
0228     // The following construct reserves stackspace for objects
0229     // of type IdDistPair_t WITHOUT initializing those objects
0230     using IdDistPair_t = HybridManager2::BoxIdDistancePair_t;
0231     char unitstackspace[VECGEOM_MAXDAUGHTERS * sizeof(IdDistPair_t)];
0232     IdDistPair_t *boxsafetylist = reinterpret_cast<IdDistPair_t *>(&unitstackspace);
0233
0234     double safety    = outsafety; // we use the outsafety estimate as starting point
0235     double safetysqr = safety * safety;
0236
0237     // safety to bounding boxes
0238     if (safety > 0. && lvol->GetDaughtersp()->size() > 0) {
0239       // calculate squared bounding box safeties in vectorized way
0240       auto ncandidates = GetSafetyCandidates_v(*fAccelerationStructureManager.GetAccStructure(lvol), localpoint,
0241                                                boxsafetylist, safetysqr);
0242       // not sorting the candidate list ( which one could do )
0243       for (unsigned int candidate = 0; candidate < ncandidates; ++candidate) {
0244         auto boxsafetypair = boxsafetylist[candidate];
0245         if (boxsafetypair.second < safetysqr) {
0246           VPlacedVolume const *cand = LookupDaughter(lvol, boxsafetypair.first);
0247           if (size_t(boxsafetypair.first) > lvol->GetDaughtersp()->size()) break;
0248           auto candidatesafety = cand->SafetyToIn(localpoint);
0249 #ifdef VERBOSE
0250           if (candidatesafety * candidatesafety > boxsafetypair.second && boxsafetypair.second > 0)
0251             std::cerr << "real safety smaller than boxsafety \n";
0252 #endif
0253           if (candidatesafety < safety) {
0254             safety    = candidatesafety;
0255             safetysqr = safety * safety;
0256           }
0257         }
0258       }
0259     }
0260     return safety;
0261   }
0262
0263   // this is (almost) the same code as in SimpleABBoxSafetyEstimator --> avoid this
0264   VECGEOM_FORCE_INLINE
0265   VECCORE_ATT_HOST_DEVICE
0266   virtual Precision ComputeSafetyForLocalPoint(Vector3D<Precision> const &localpoint,
0267                                                VPlacedVolume const *pvol) const override
0268   {
0269     // safety to mother
0270     double safety = pvol->SafetyToOut(localpoint);
0271     if (safety <= 0.) {
0272       return 0.;
0273     }
0274     return TreatSafetyToIn(localpoint, pvol->GetLogicalVolume(), safety);
0275   }
0276
0277   VECGEOM_FORCE_INLINE
0278   VECCORE_ATT_HOST_DEVICE
0279   virtual Precision ComputeSafetyToDaughtersForLocalPoint(Vector3D<Precision> const &localpoint,
0280                                                           LogicalVolume const *lvol) const override
0281   {
0282     return TreatSafetyToIn(localpoint, lvol, kInfLength);
0283   }
0284
0285   // vector interface
0286   VECGEOM_FORCE_INLINE
0287   virtual void ComputeSafetyForLocalPoints(SOA3D<Precision> const & /*localpoints*/, VPlacedVolume const * /*pvol*/,
0288                                            Precision * /*safeties*/) const override
0289   {
0290     //    // a stack based workspace array
0291     //    static __thread ABBoxManager::BoxIdDistancePair_t boxsafetylist[VECGEOM_MAXDAUGHTERS] = {};
0292     //
0293     //    // safety to mother -- using vector interface
0294     //    pvol->SafetyToOut(localpoints, safeties);
0295     //
0296     //    // safety to bounding boxes
0297     //    LogicalVolume const *lvol = pvol->GetLogicalVolume();
0298     //    if (!(lvol->GetDaughtersp()->size() > 0))
0299     //      return;
0300     //
0301     //    // get bounding boxes (they are the same for all tracks)
0302     //    int numberofboxes;
0303     //    auto bboxes = fABBoxManager.GetABBoxes_v(lvol, numberofboxes);
0304     //
0305     //    // now loop over particles
0306     //    for (int i = 0, ntracks = localpoints.size(); i < ntracks; ++i) {
0307     //      double safety = safeties[i];
0308     //      if (safeties[i] > 0.) {
0309     //        double safetysqr = safeties[i] * safeties[i];
0310     //        auto lpoint = localpoints[i];
0311     //        // vectorized search through bounding boxes -- quickly excluding many candidates
0312     //        auto ncandidates = GetSafetyCandidates_v(lpoint, bboxes, numberofboxes, boxsafetylist, safetysqr);
0313     //        // loop over remaining candidates
0314     //        for (unsigned int candidate = 0; candidate < ncandidates; ++candidate) {
0315     //          auto boxsafetypair = boxsafetylist[candidate];
0316     //          if (boxsafetypair.second < safetysqr) {
0317     //            VPlacedVolume const *candidate = LookupDaughter(lvol, boxsafetypair.first);
0318     //            if (boxsafetypair.first > lvol->GetDaughtersp()->size())
0319     //              break;
0320     //            auto candidatesafety = candidate->SafetyToIn(lpoint);
0321     //#ifdef VERBOSE
0322     //            if (candidatesafety * candidatesafety > boxsafetypair.second && boxsafetypair.second > 0)
0323     //              std::cerr << "real safety smaller than boxsafety \n";
0324     //#endif
0325     //            if (candidatesafety < safety) {
0326     //              safety = candidatesafety;
0327     //              safetysqr = safety * safety;
0328     //            }
0329     //          }
0330     //        }
0331     //      }
0332     //      // write back result
0333     //      safeties[i] = safety;
0334     //    }
0335   }
0336
0337   static VSafetyEstimator *Instance()
0338   {
0339     static HybridSafetyEstimator instance;
0340     return &instance;
0341   }
0342
0343 }; // end class
0344 } // namespace VECGEOM_IMPL_NAMESPACE
0345 } // namespace vecgeom
0346
0347 #endif /* NAVIGATION_SIMPLEABBOXSAFETYESTIMATOR_H_ */