Commit 2cc22f2c authored by Sandro Christian Wenzel's avatar Sandro Christian Wenzel
Browse files

Prototypic versions of SafetyEstimator and LevelLocator using voxel structure

Introducing navigation classes that estimate safety and do the "level locate"
operation using the voxel structure managed by FlatVoxelManager.

In the large daughter limit, this provides better scaling than just using
the estimators based on BVH structure, since "safety" and "locate"
are very local operations.

This is a very first version for demonstration purposes.

Future directions/work will involve:
  a) automatic tuning of the right voxel size/structure
  b) combination of voxel structure with SIMD BVH
  c) speeding up creation of voxel structure
     (it is multi-threaded but we can make use of SIMD and GPU)

Adapt benchmark for SafetyEstimator accordingly.
parent 9568e9bd
// This file is part of VecGeom and is distributed under the
// conditions in the file LICENSE.txt in the top directory.
// For the full list of authors see CONTRIBUTORS.txt and `git log`.
/// \file navigation/VoxelLevelLocator.h
/// \author Sandro Wenzel (CERN)
#ifndef NAVIGATION_VOXELLEVELLOCATOR_H_
#define NAVIGATION_VOXELLEVELLOCATOR_H_
#include "navigation/VLevelLocator.h"
#include "volumes/LogicalVolume.h"
#include "volumes/PlacedVolume.h"
#include "management/FlatVoxelManager.h"
#include "volumes/kernel/BoxImplementation.h"
#include "navigation/SimpleLevelLocator.h"
#include "management/ABBoxManager.h"
namespace vecgeom {
inline namespace VECGEOM_IMPL_NAMESPACE {
//! A LevelLocator using voxel hash-maps containing candidate volume ids to quickly
//! decide in which volume a point is.
template <bool IsAssemblyAware = false>
class TVoxelLevelLocator : public VLevelLocator {
private:
FlatVoxelManager &fAccelerationStructure;
TVoxelLevelLocator() : fAccelerationStructure(FlatVoxelManager::Instance()) {}
// the actual implementation kernel
// the template "ifs" should be optimized away
// arguments are pointers to allow for nullptr
template <bool ExclV, bool ModifyState>
__attribute__((always_inline)) bool LevelLocateKernel(LogicalVolume const *lvol, VPlacedVolume const *exclvol,
Vector3D<Precision> const &localpoint, NavigationState *state,
VPlacedVolume const *&pvol,
Vector3D<Precision> &daughterlocalpoint) const
{
const auto accstructure = fAccelerationStructure.GetStructure(lvol);
// fetch the right voxels
const auto locatevoxels = accstructure->fVoxelToLocateCandidates;
// fetch the list of candidates to check (as id)
Vector3D<float> floatlocalpoint(localpoint.x(), localpoint.y(), localpoint.z());
int numbercandidates{0};
const auto candidateidsptr = locatevoxels->getProperties(floatlocalpoint, numbercandidates);
if (numbercandidates > 0) {
// here we have something to do
int numberboxes{0};
const auto boxes = ABBoxManager::Instance().GetABBoxes(lvol, numberboxes);
for (int i = 0; i < numbercandidates; ++i) {
const int daughterid = candidateidsptr[i];
// discard blocked volume
VPlacedVolume const *candidate = lvol->GetDaughters()[daughterid];
if (ExclV) {
if (candidate == exclvol) continue;
}
// quick aligned bounding box check (could be done in SIMD fashion if multiple candidates)
const auto &lower = boxes[2 * daughterid];
const auto &upper = boxes[2 * daughterid + 1];
bool boxcontains{false};
ABBoxImplementation::ABBoxContainsKernel(lower, upper, localpoint, boxcontains);
if (!boxcontains) {
continue;
}
// final candidate check
if (CheckCandidateVol<IsAssemblyAware, ModifyState>(candidate, localpoint, state, pvol, daughterlocalpoint)) {
return true;
}
}
}
return false;
}
// the actual implementation kernel
// the template "ifs" should be optimized away
// arguments are pointers to allow for nullptr
template <bool ExclV, bool ModifyState>
__attribute__((always_inline)) bool LevelLocateKernelWithDirection(LogicalVolume const *lvol,
VPlacedVolume const *exclvol,
Vector3D<Precision> const &localpoint,
Vector3D<Precision> const &localdir,
NavigationState *state, VPlacedVolume const *&pvol,
Vector3D<Precision> &daughterlocalpoint) const
{
const auto accstructure = fAccelerationStructure.GetStructure(lvol);
// fetch the right voxels
const auto locatevoxels = accstructure->fVoxelToLocateCandidates;
// fetch the list of candidates to check (as id)
Vector3D<float> floatlocalpoint(localpoint.x(), localpoint.y(), localpoint.z());
int numbercandidates{0};
const auto candidateidsptr = locatevoxels->getProperties(floatlocalpoint, numbercandidates);
if (numbercandidates > 0) {
// here we have something to do
int numberboxes{0};
const auto boxes = ABBoxManager::Instance().GetABBoxes(lvol, numberboxes);
for (int i = 0; i < numbercandidates; ++i) {
const int daughterid = candidateidsptr[i];
// discard blocked volume
VPlacedVolume const *candidate = lvol->GetDaughters()[daughterid];
if (ExclV) {
if (candidate == exclvol) continue;
}
// quick aligned bounding box check (could be done in SIMD fashion if multiple candidates)
const auto &lower = boxes[2 * daughterid];
const auto &upper = boxes[2 * daughterid + 1];
bool boxcontains{false};
ABBoxImplementation::ABBoxContainsKernel(lower, upper, localpoint, boxcontains);
if (!boxcontains) {
continue;
}
// final candidate check
if (CheckCandidateVolWithDirection<IsAssemblyAware, ModifyState>(candidate, localpoint, localdir, state, pvol,
daughterlocalpoint)) {
return true;
}
}
}
return false;
}
public:
static std::string GetClassName() { return "VoxelLevelLocator"; }
std::string GetName() const override { return GetClassName(); }
bool LevelLocate(LogicalVolume const *lvol, Vector3D<Precision> const &localpoint, VPlacedVolume const *&pvol,
Vector3D<Precision> &daughterlocalpoint) const override
{
return LevelLocateKernel<false, false>(lvol, nullptr, localpoint, nullptr, pvol, daughterlocalpoint);
}
bool LevelLocate(LogicalVolume const *lvol, Vector3D<Precision> const &localpoint, NavigationState &state,
Vector3D<Precision> &daughterlocalpoint) const override
{
VPlacedVolume const *pvol;
return LevelLocateKernel<false, true>(lvol, nullptr, localpoint, &state, pvol, daughterlocalpoint);
}
bool LevelLocateExclVol(LogicalVolume const *lvol, VPlacedVolume const *exclvol,
Vector3D<Precision> const &localpoint, VPlacedVolume const *&pvol,
Vector3D<Precision> &daughterlocalpoint) const override
{
return LevelLocateKernel<true, false>(lvol, exclvol, localpoint, nullptr, pvol, daughterlocalpoint);
}
bool LevelLocateExclVol(LogicalVolume const *lvol, VPlacedVolume const *exclvol,
Vector3D<Precision> const &localpoint, Vector3D<Precision> const &localdir,
VPlacedVolume const *&pvol, Vector3D<Precision> &daughterlocalpoint) const override
{
return LevelLocateKernelWithDirection<true, false>(lvol, exclvol, localpoint, localdir, nullptr, pvol,
daughterlocalpoint);
}
static VLevelLocator const *GetInstance()
{
static TVoxelLevelLocator instance;
return &instance;
}
}; // end class declaration
template <>
inline std::string TVoxelLevelLocator<true>::GetClassName()
{
return "AssemblyAwareVoxelLevelLocator";
}
} // namespace VECGEOM_IMPL_NAMESPACE
} // namespace vecgeom
#endif /* NAVIGATION_VOXELLEVELLOCATOR_H_ */
// This file is part of VecGeom and is distributed under the
// conditions in the file LICENSE.txt in the top directory.
// For the full list of authors see CONTRIBUTORS.txt and `git log`.
/// \file navigation/VoxelSafetyEstimator.h
/// \author Sandro Wenzel (CERN)
#ifndef NAVIGATION_VOXELSAFETYESTIMATOR_H_
#define NAVIGATION_VOXELSAFETYESTIMATOR_H_
#include "navigation/VSafetyEstimator.h"
#include "management/FlatVoxelManager.h"
#include <cmath>
namespace vecgeom {
inline namespace VECGEOM_IMPL_NAMESPACE {
//! a safety estimator using a voxel lookup table of precomputed safety candidates
//! (or in the extreme case of precomputed safety values)
class VoxelSafetyEstimator : public VSafetyEstimatorHelper<VoxelSafetyEstimator> {
private:
const FlatVoxelManager &fAccStructureManager;
VoxelSafetyEstimator()
: VSafetyEstimatorHelper<VoxelSafetyEstimator>(), fAccStructureManager{FlatVoxelManager::Instance()}
{
}
// convert index to physical daugher
VPlacedVolume const *LookupDaughter(LogicalVolume const *lvol, int id) const
{
assert(id >= 0 && "access with negative index");
assert(size_t(id) < lvol->GetDaughtersp()->size() && "access beyond size of daughterlist");
return lvol->GetDaughtersp()->operator[](id);
}
public:
static constexpr const char *gClassNameString = "VoxelSafetyEstimator";
VECGEOM_FORCE_INLINE
VECCORE_ATT_HOST_DEVICE
Precision TreatSafetyToIn(Vector3D<Precision> const &localpoint, LogicalVolume const *lvol, Precision outsafety) const
{
throw std::runtime_error("unimplemented function called");
}
// #define CROSSCHECK
#ifdef CROSSCHECK
void printProcedure(int const *safetycandidates, int length, LogicalVolume const *lvol,
Vector3D<Precision> const &localpoint) const
{
int size{0};
auto abboxcorners = ABBoxManager::Instance().GetABBoxes(lvol, size);
Vector3D<float> lp(localpoint.x(), localpoint.y(), localpoint.z()); // in float
const bool needmother = (safetycandidates[0] == -1);
const Precision safetymother = needmother ? lvol->GetUnplacedVolume()->SafetyToOut(localpoint) : 0.;
std::cerr << "safetymother " << safetymother;
Precision safetyToDaughters{1E20};
// now talk to daughter candidates - with bounding box speedup
// TODO: this could actually happen with SIMD speedup
float bestSafetySqr{1E20};
for (int i = 0; i < length; ++i) {
if (safetycandidates[i] == -1) continue;
const auto candidateid = safetycandidates[i];
const auto &lower = abboxcorners[2 * candidateid];
const auto &upper = abboxcorners[2 * candidateid + 1];
const float ssqr = ABBoxImplementation::ABBoxSafetySqr(lower, upper, lp);
std::cerr << i << " boxsqr " << ssqr << "\t" << std::sqrt(ssqr) << "\n";
if (ssqr > 0) {
bestSafetySqr = std::min(bestSafetySqr, ssqr);
} else { // inside
const auto daughter = LookupDaughter(lvol, safetycandidates[i]);
auto sin = daughter->SafetyToIn(localpoint);
std::cerr << i << " SI " << sin << "\n";
safetyToDaughters = std::min(safetyToDaughters, daughter->SafetyToIn(localpoint));
}
}
safetyToDaughters = std::min(safetyToDaughters, (Precision)std::sqrt(bestSafetySqr));
std::cerr << "final safetyd " << safetyToDaughters << "\n";
const float returnvalue = needmother ? std::min(safetyToDaughters, safetymother) : safetyToDaughters;
std::cerr << "returning " << returnvalue << "\n";
}
#endif
VECGEOM_FORCE_INLINE
VECCORE_ATT_HOST_DEVICE
virtual Precision ComputeSafetyForLocalPoint(Vector3D<Precision> const &localpoint,
VPlacedVolume const *pvol) const override
{
static int counter = 0;
counter++;
const auto lvol = pvol->GetLogicalVolume();
const auto structure = fAccStructureManager.GetStructure(lvol);
int size{0};
auto abboxcorners = ABBoxManager::Instance().GetABBoxes(lvol, size);
if (structure != nullptr) {
const Vector3D<float> lp(localpoint.x(), localpoint.y(), localpoint.z()); // in float
// fetch safety properties for localpoint;
const int *safetycandidates{nullptr};
int length{0};
safetycandidates = structure->fVoxelToCandidate->getProperties(lp, length);
if (length > 0) {
const bool needmother = true; //(safetycandidates[0] == -1);
const Precision safetymother = needmother ? lvol->GetUnplacedVolume()->SafetyToOut(localpoint) : 0.;
Precision safetyToDaughters{1E20};
// now talk to daughter candidates - with bounding box speedup
// TODO: this could actually happen with SIMD speedup
float bestSafetySqr{1E20};
for (int i = 0; i < length; ++i) {
if (safetycandidates[i] == -1) continue;
const auto candidateid = safetycandidates[i];
const auto &lower = abboxcorners[2 * candidateid];
const auto &upper = abboxcorners[2 * candidateid + 1];
const float ssqr = ABBoxImplementation::ABBoxSafetySqr(lower, upper, lp);
if (ssqr > 0) {
bestSafetySqr = std::min(bestSafetySqr, ssqr);
} else { // inside
const auto daughter = LookupDaughter(lvol, safetycandidates[i]);
auto s = daughter->SafetyToIn(localpoint);
if (s <= 0.) {
// can break here
s = 0.;
}
safetyToDaughters = std::min(safetyToDaughters, s);
}
}
safetyToDaughters = std::min(safetyToDaughters, (Precision)std::sqrt(bestSafetySqr));
const float returnvalue = needmother ? std::min(safetyToDaughters, safetymother) : safetyToDaughters;
#ifdef CROSSCHECK
if (returnvalue > lvol->GetUnplacedVolume()->SafetyToOut(localpoint) + 1E-4) {
std::cerr << returnvalue << " " << lvol->GetUnplacedVolume()->SafetyToOut(localpoint)
<< "Problem in voxel safety for point " << lp << " voxelkey "
<< structure->fVoxelToCandidate->getKey(lp.x(), lp.y(), lp.z()) << " ";
std::cerr << "{ ";
for (int i = 0; i < length; ++i) {
std::cout << safetycandidates[i] << " , ";
}
std::cerr << " }\n";
printProcedure(safetycandidates, length, lvol, localpoint);
}
#endif
if (returnvalue < 0) {
std::cerr << "returning negative value " << safetyToDaughters << " " << safetymother << "\n";
}
return returnvalue;
} else {
// no information for this voxel present
std::cerr << counter << " no information for this voxel present " << localpoint << " at key "
<< structure->fVoxelToCandidate->getKey(lp.x(), lp.y(), lp.z()) << " \n ";
return 0.;
}
}
std::cerr << " No voxel information found; Cannot use voxel safety\n";
return 0.;
}
VECGEOM_FORCE_INLINE
VECCORE_ATT_HOST_DEVICE
virtual Precision ComputeSafetyToDaughtersForLocalPoint(Vector3D<Precision> const &localpoint,
LogicalVolume const *lvol) const override
{
return TreatSafetyToIn(localpoint, lvol, kInfLength);
}
//
// These interfaces for baskets are only dummy-implemented at the moment
//
VECCORE_ATT_HOST_DEVICE
Real_v ComputeSafetyForLocalPoint(Vector3D<Real_v> const & /*localpoint*/, VPlacedVolume const * /*pvol*/,
Bool_v /*m*/) const override
{
throw std::runtime_error("unimplemented function called");
return Real_v(0.);
}
// interfaces to treat vectors/collections of points (uses the approach with intermediate storage and passing down the
// loops to shapes)
void ComputeVectorSafety(SOA3D<Precision> const & /*globalpoints*/, NavStatePool &states,
SOA3D<Precision> & /*workspace*/, Precision * /*safeties*/) const override
{
throw std::runtime_error("unimplemented function called");
}
// interfaces to treat vectors/collections of points (uses the approach without intermediate storage; requires access
// to new SIMD interface)
void ComputeVectorSafety(SOA3D<Precision> const & /*globalpoints*/, NavStatePool & /*states*/,
Precision * /*safeties*/) const override
{
throw std::runtime_error("unimplemented function called");
}
void ComputeSafetyForLocalPoints(SOA3D<Precision> const & /*localpoints*/, VPlacedVolume const * /*pvol*/,
Precision * /*safeties*/) const override
{
throw std::runtime_error("unimplemented function called");
}
static VSafetyEstimator *Instance()
{
static VoxelSafetyEstimator instance;
return &instance;
}
}; // end class
} // namespace VECGEOM_IMPL_NAMESPACE
} // namespace vecgeom
#endif /* NAVIGATION_VOXELSAFETYESTIMATOR_H_ */
......@@ -408,8 +408,8 @@ FlatVoxelHashMap<int, false> *FlatVoxelManager::BuildSafetyVoxels(LogicalVolume
}
std::cerr << " done \n";
safetyvoxels->print();
// create cash
// safetyvoxels->print();
// create cache
safetyvoxels->dumpToTFile(createName(vol, Nx, Ny, Nz).c_str());
return safetyvoxels;
......@@ -495,7 +495,7 @@ FlatVoxelHashMap<int, false> *FlatVoxelManager::BuildLocateVoxels(LogicalVolume
locatevoxels->addPropertyForKey(key, cand);
}
}
locatevoxels->print();
// locatevoxels->print();
return locatevoxels;
}
......
......@@ -22,6 +22,8 @@
#include "navigation/SimpleABBoxSafetyEstimator.h"
#include "navigation/HybridSafetyEstimator.h"
#include "management/HybridManager2.h"
#include "navigation/VoxelSafetyEstimator.h"
#include "management/FlatVoxelManager.h"
// in case someone has written a special safety estimator for the CAHL logical volume
//#include "navigation/CAHLSafetyEstimator.h"
......@@ -35,6 +37,14 @@
#include "TGeoBBox.h"
#endif
#ifdef VECGEOM_GEANT4
#include "G4Navigator.hh"
#include "G4VPhysicalVolume.hh"
#include "G4ThreeVector.hh"
#include "management/G4GeoManager.h"
#include "G4VoxelNavigation.hh"
#endif
#include <iostream>
//#define CALLGRIND_ENABLED
......@@ -75,6 +85,77 @@ __attribute__((noinline)) void benchSafety(SOA3D<Precision> const &points, NavSt
std::cerr << "accum " << T::GetClassName() << " " << accum << "\n";
}
template <typename T>
__attribute__((noinline)) void benchLocalSafety(SOA3D<Precision> const &localpoints, NavStatePool &pool)
{
// bench safety
Precision *safety = new Precision[localpoints.size()];
Stopwatch timer;
VSafetyEstimator *se = T::Instance();
const auto topvolume = pool[0]->Top();
timer.Start();
for (size_t i = 0; i < localpoints.size(); ++i) {
safety[i] = se->ComputeSafetyForLocalPoint(localpoints[i], topvolume);
}
timer.Stop();
std::cerr << timer.Elapsed() << "\n";
double accum(0.);
for (size_t i = 0; i < localpoints.size(); ++i) {
accum += safety[i];
}
delete[] safety;
std::cerr << "accum " << T::GetClassName() << " " << accum << "\n";
}
#ifdef VECGEOM_GEANT4
__attribute__((noinline)) void benchmarkLocalG4Safety(SOA3D<Precision> const &points,
SOA3D<Precision> const &localpoints)
{
Stopwatch timer;
G4VPhysicalVolume *world(vecgeom::G4GeoManager::Instance().GetG4GeometryFromROOT());
if (world != nullptr) G4GeoManager::Instance().LoadG4Geometry(world);
// Note: Vector3D's are expressed in cm, while G4ThreeVectors are expressed in mm
const Precision cm = 10.; // cm --> mm conversion
G4Navigator &g4nav = *(G4GeoManager::Instance().GetNavigator());
Precision *safety = new Precision[points.size()];
G4VoxelNavigation g4voxelnavigator;
// we need one navigation history object for the points; taking the first point is enough
// as all of them are in the same toplevel volume
G4ThreeVector g4pos(points[0].x() * cm, points[0].y() * cm, points[0].z() * cm);
g4nav.LocateGlobalPointAndSetup(g4pos, nullptr, false);
auto history = g4nav.CreateTouchableHistory()->GetHistory();
std::cerr << history << "\n";
std::cerr << history->GetDepth() << "\n";
#ifdef CALLGRIND_ENABLED
CALLGRIND_START_INSTRUMENTATION;
#endif
timer.Start();
for (decltype(points.size()) i = 0; i < points.size(); ++i) {
const auto &vgpoint = points[i];
const G4ThreeVector g4lpos(vgpoint.x() * cm, vgpoint.y() * cm, vgpoint.z() * cm);
safety[i] = g4nav.ComputeSafety(g4lpos); // *history, 1E20);
}
timer.Stop();
std::cerr << (Precision)timer.Elapsed() << "\n";
#ifdef CALLGRIND_ENABLED
CALLGRIND_STOP_INSTRUMENTATION;
CALLGRIND_DUMP_STATS;
#endif
double accum(0.);
for (size_t i = 0; i < localpoints.size(); ++i) {
accum += safety[i];
}
// cleanup
delete[] safety;
std::cerr << "accum G4 " << accum / cm << "\n";
}
#endif // end if G4
#ifdef VECGEOM_VECTOR
// benchmarks the old vector interface (needing temporary workspace)
template <typename T>
__attribute__((noinline)) void benchVectorSafety(SOA3D<Precision> const &points, NavStatePool &pool)
......@@ -115,12 +196,12 @@ __attribute__((noinline)) void benchVectorSafetyNoWorkspace(SOA3D<Precision> con
vecCore::AlignedFree(safety);
std::cerr << "VECTOR (NO WORKSP) accum " << T::GetClassName() << " " << accum << "\n";
}
#endif
// benchmarks the ROOT safety interface
#ifdef VECGEOM_ROOT
__attribute__((noinline)) void benchmarkROOTSafety(int nPoints, SOA3D<Precision> const &points, int i)
__attribute__((noinline)) void benchmarkROOTSafety(int nPoints, SOA3D<Precision> const &points)
{
TGeoNavigator *rootnav = ::gGeoManager->GetCurrentNavigator();
TGeoBranchArray *brancharrays[nPoints];
Precision *safety = new Precision[nPoints];
......@@ -155,32 +236,90 @@ __attribute__((noinline)) void benchmarkROOTSafety(int nPoints, SOA3D<Precision>
accum += safety[i];
}
std::cerr << "ROOT s " << accum << "\n";
delete[] safety;
return;
}
__attribute__((noinline)) void benchmarkLocalROOTSafety(int nPoints, SOA3D<Precision> const &points,
SOA3D<Precision> const &localpoints)
{
// we init the ROOT navigator to a specific branch state
auto nav = gGeoManager->GetCurrentNavigator();
nav->FindNode(points[0].x(), points[0].y(), points[0].z());
Precision *safety = new Precision[nPoints];
#ifdef CALLGRIND_ENABLED
CALLGRIND_START_INSTRUMENTATION;
#endif
Stopwatch timer;
timer.Start();
for (int i = 0; i < nPoints; ++i) {
// There is no direct way of calling the local safety function; points are always transformed so I have to
// give the global point
Vector3D<Precision> const &pos = points[i];
nav->SetCurrentPoint(pos.x(), pos.y(), pos.z());
safety[i] = nav->Safety();
}
timer.Stop();
std::cerr << "ROOT time" << timer.Elapsed() << "\n";
#ifdef CALLGRIND_ENABLED
CALLGRIND_STOP_INSTRUMENTATION;
CALLGRIND_DUMP_STATS;
#endif
double accum(0.);
for (int i = 0; i < nPoints; ++i) {
accum += safety[i];
}
std::cerr << "ROOT s " << accum << "\n";
delete[] safety;
return;
}
#endif
// main routine starting up the individual benchmarks
void benchDifferentSafeties(SOA3D<Precision> const &points, NavStatePool &pool)
void benchDifferentSafeties(SOA3D<Precision> const &points, SOA3D<Precision> const &localpoints, NavStatePool &pool)
{
std::cerr << "##\n";
std::cerr << "## - GLOBAL POINTS - \n";
std::cerr << "##\n";
RUNBENCH(benchSafety<SimpleSafetyEstimator>(points, pool));
std::cerr << "##\n";
RUNBENCH(benchSafety<SimpleABBoxSafetyEstimator>(points, pool));
std::cerr << "##\n";
RUNBENCH(benchSafety<HybridSafetyEstimator>(points, pool));
std::cerr << "##\n";
RUNBENCH(benchSafety<VoxelSafetyEstimator>(points, pool));