diff --git a/.gitignore b/.gitignore index c0a676728ea..4afa640a655 100644 --- a/.gitignore +++ b/.gitignore @@ -35,3 +35,9 @@ build* !.gitlab-ci.yml !.clang-tidy !.oclint + +Examples/Scripts/Python/csv +Examples/Scripts/Python/onnx_models +Examples/Scripts/Python/torchscript_models/ +*.old +*.old2 diff --git a/CMakeLists.txt b/CMakeLists.txt index 946452b2bd3..be9102cdb62 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -42,6 +42,8 @@ option(ACTS_BUILD_PLUGIN_GEANT4 "Build Geant4 plugin" OFF) option(ACTS_BUILD_PLUGIN_EXATRKX "Build the Exa.TrkX plugin" OFF) option(ACTS_EXATRKX_ENABLE_ONNX "Build the Onnx backend for the exatrkx plugin" OFF) option(ACTS_EXATRKX_ENABLE_TORCH "Build the torchscript backend for the exatrkx plugin" ON) +option(ACTS_EXATRKX_ENABLE_STUB "Only build a stub backend for the exatrkx plugin" OFF) +option(ACTS_USE_SYSTEM_ACTSDD4HEP "Use the ActsDD4hep glue library provided by the system instead of building it" OFF) option(ACTS_BUILD_PLUGIN_IDENTIFICATION "Build Identification plugin" OFF) option(ACTS_BUILD_PLUGIN_JSON "Build json plugin" OFF) option(ACTS_USE_SYSTEM_NLOHMANN_JSON "Use nlohmann::json provided by the system instead of the bundled version" ${ACTS_USE_SYSTEM_LIBS}) diff --git a/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp b/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp index 930f86dd7a5..95e4a708152 100644 --- a/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp +++ b/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp @@ -312,7 +312,7 @@ class CombinatorialKalmanFilter { m_updaterLogger{m_logger->cloneWithSuffix("Updater")}, m_smootherLogger{m_logger->cloneWithSuffix("Smoother")} {} - private: + public: using KalmanNavigator = typename propagator_t::Navigator; /// The propagator for the transport and material update diff --git a/Core/include/Acts/Utilities/ContainerPrinter.hpp b/Core/include/Acts/Utilities/ContainerPrinter.hpp new file mode 100644 index 00000000000..a28ce84299f --- /dev/null +++ b/Core/include/Acts/Utilities/ContainerPrinter.hpp @@ -0,0 +1,53 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2023 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +#include +#include + +namespace Acts { + +template +struct ContainerPrinter { + iterator_t begin; + iterator_t end; + + template + ContainerPrinter(const container_t &c, std::size_t max) { + begin = c.cbegin(); + const auto n = std::min(c.size(), max); + end = begin; + std::advance(end, n); + assert(std::distance(begin, end) == static_cast(n)); + } + + template + ContainerPrinter(const container_t &c) : ContainerPrinter(c, c.size()) {} + + ContainerPrinter(iterator_t a, iterator_t b) : begin(a), end(b) {} +}; + +template +std::ostream &operator<<(std::ostream &os, + const ContainerPrinter &p) { + for (auto it = p.begin; it != p.end; ++it) { + os << *it << " "; + } + return os; +} + +template +ContainerPrinter(const container_t &, std::size_t) + -> ContainerPrinter::const_iterator>; + +template +ContainerPrinter(const container_t &) + -> ContainerPrinter::const_iterator>; + +} // namespace Acts diff --git a/Core/include/Acts/Utilities/KDTree.hpp b/Core/include/Acts/Utilities/KDTree.hpp index d78f8cf021c..0cae5c87d27 100644 --- a/Core/include/Acts/Utilities/KDTree.hpp +++ b/Core/include/Acts/Utilities/KDTree.hpp @@ -12,6 +12,7 @@ #include #include +#include #include #include #include diff --git a/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationAlgorithm.hpp b/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationAlgorithm.hpp index 443390583e7..5d74efe6760 100644 --- a/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationAlgorithm.hpp +++ b/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationAlgorithm.hpp @@ -26,8 +26,7 @@ #include "ActsExamples/Framework/ProcessCode.hpp" #include "ActsExamples/Framework/RandomNumbers.hpp" #include "ActsFatras/Digitization/Channelizer.hpp" -#include "ActsFatras/Digitization/PlanarSurfaceDrift.hpp" -#include "ActsFatras/Digitization/PlanarSurfaceMask.hpp" +#include "ActsFatras/Digitization/Segmentizer.hpp" #include "ActsFatras/Digitization/UncorrelatedHitSmearer.hpp" #include @@ -69,20 +68,6 @@ class DigitizationAlgorithm final : public IAlgorithm { const DigitizationConfig& config() const { return m_cfg; } private: - /// Helper method for the geometric channelizing part - /// - /// @param geoCfg is the geometric digitization configuration - /// @param hit the Simultated hit - /// @param surface the Surface on which this is supposed to happen - /// @param gctx the Geometry context - /// @param rng the Random number engine for the drift smearing - /// - /// @return the list of channels - std::vector channelizing( - const GeometricConfig& geoCfg, const SimHit& hit, - const Acts::Surface& surface, const Acts::GeometryContext& gctx, - RandomEngine& rng) const; - /// Helper method for creating digitized parameters from clusters /// /// @todo ADD random smearing @@ -93,7 +78,7 @@ class DigitizationAlgorithm final : public IAlgorithm { /// @return the list of digitized parameters DigitizedParameters localParameters( const GeometricConfig& geoCfg, - const std::vector& channels, + const std::vector& channels, RandomEngine& rng) const; /// Nested smearer struct that holds geometric digitizer and smearing @@ -113,9 +98,7 @@ class DigitizationAlgorithm final : public IAlgorithm { DigitizationConfig m_cfg; /// Digitizers within geometry hierarchy Acts::GeometryHierarchyMap m_digitizers; - /// Geometric digtizers - ActsFatras::PlanarSurfaceDrift m_surfaceDrift; - ActsFatras::PlanarSurfaceMask m_surfaceMask; + /// Geometric digtizer ActsFatras::Channelizer m_channelizer; ReadDataHandle m_simContainerReadHandle{this, diff --git a/Examples/Algorithms/Digitization/src/DigitizationAlgorithm.cpp b/Examples/Algorithms/Digitization/src/DigitizationAlgorithm.cpp index cc2fa9bb790..dc3df736766 100644 --- a/Examples/Algorithms/Digitization/src/DigitizationAlgorithm.cpp +++ b/Examples/Algorithms/Digitization/src/DigitizationAlgorithm.cpp @@ -205,16 +205,20 @@ ActsExamples::ProcessCode ActsExamples::DigitizationAlgorithm::execute( ACTS_VERBOSE("Configured to geometric digitize " << digitizer.geometric.indices.size() << " parameters."); - auto channels = channelizing(digitizer.geometric, simHit, - *surfacePtr, ctx.geoContext, rng); - if (channels.empty()) { + const auto& cfg = digitizer.geometric; + Acts::Vector3 driftDir = cfg.drift(simHit.position(), rng); + auto channelsRes = m_channelizer.channelize( + simHit, *surfacePtr, ctx.geoContext, driftDir, + cfg.segmentation, cfg.thickness); + if (!channelsRes.ok() || channelsRes->empty()) { ACTS_DEBUG( "Geometric channelization did not work, skipping this hit.") continue; } - ACTS_VERBOSE("Activated " << channels.size() + ACTS_VERBOSE("Activated " << channelsRes->size() << " channels for this hit."); - dParameters = localParameters(digitizer.geometric, channels, rng); + dParameters = + localParameters(digitizer.geometric, *channelsRes, rng); } // Smearing part - (optionally) rest @@ -292,30 +296,10 @@ ActsExamples::ProcessCode ActsExamples::DigitizationAlgorithm::execute( return ProcessCode::SUCCESS; } -std::vector -ActsExamples::DigitizationAlgorithm::channelizing( - const GeometricConfig& geoCfg, const SimHit& hit, - const Acts::Surface& surface, const Acts::GeometryContext& gctx, - RandomEngine& rng) const { - Acts::Vector3 driftDir = geoCfg.drift(hit.position(), rng); - - auto driftedSegment = - m_surfaceDrift.toReadout(gctx, surface, geoCfg.thickness, hit.position(), - hit.direction(), driftDir); - auto maskedSegmentRes = m_surfaceMask.apply(surface, driftedSegment); - if (maskedSegmentRes.ok()) { - auto maskedSegment = maskedSegmentRes.value(); - // Now Channelize - return m_channelizer.segments(gctx, surface, geoCfg.segmentation, - maskedSegment); - } - return {}; -} - ActsExamples::DigitizedParameters ActsExamples::DigitizationAlgorithm::localParameters( const GeometricConfig& geoCfg, - const std::vector& channels, + const std::vector& channels, RandomEngine& rng) const { DigitizedParameters dParameters; diff --git a/Examples/Algorithms/Geant4/src/SensitiveSteppingAction.cpp b/Examples/Algorithms/Geant4/src/SensitiveSteppingAction.cpp index a574695c1eb..30a7e5ea9ff 100644 --- a/Examples/Algorithms/Geant4/src/SensitiveSteppingAction.cpp +++ b/Examples/Algorithms/Geant4/src/SensitiveSteppingAction.cpp @@ -28,6 +28,7 @@ #include #include #include +#include class G4PrimaryParticle; diff --git a/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackFindingAlgorithm.hpp b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackFindingAlgorithm.hpp index 4add3628212..e3279471529 100644 --- a/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackFindingAlgorithm.hpp +++ b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackFindingAlgorithm.hpp @@ -141,6 +141,10 @@ class TrackFindingAlgorithm final : public IAlgorithm { mutable std::atomic m_nTotalSeeds{0}; mutable std::atomic m_nFailedSeeds{0}; + mutable std::mutex m_mutex; + mutable std::vector m_nTracksPerSeeds; + mutable std::vector m_nSelTracksPerSeeds; + mutable tbb::combinable m_memoryStatistics{[]() { auto mtj = std::make_shared(); diff --git a/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFindingX/CombinedKfAndCkf.hpp b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFindingX/CombinedKfAndCkf.hpp new file mode 100644 index 00000000000..52718d6f60c --- /dev/null +++ b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFindingX/CombinedKfAndCkf.hpp @@ -0,0 +1,137 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2022 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/TrackFinding/CombinatorialKalmanFilter.hpp" +#include "Acts/TrackFitting/KalmanFitter.hpp" + +namespace Acts { + +template +struct CombinedKfAndCkf { + propagator_t m_propagator; + std::unique_ptr m_logger; + KalmanFitter m_kalmanFitter; + + CombinedKfAndCkf(propagator_t pPropagator, + std::unique_ptr _logger = + getDefaultLogger("KfCkfComb", Logging::INFO)) + : m_propagator(std::move(pPropagator)), + m_logger(std::move(_logger)), + m_kalmanFitter(pPropagator, m_logger->cloneWithSuffix(":KF")) {} + + const Logger& logger() const { return *m_logger; } + + template + auto runKalmanFitter( + sli_kf_t it, sli_kf_t end, const start_parameters_t& sParameters, + const CombinatorialKalmanFilterOptions& ckfOptions, + track_container_t& trackContainer) const { + KalmanFitterExtensions extensions; + extensions.calibrator = ckfOptions.extensions.calibrator; + extensions.updater = ckfOptions.extensions.updater; + extensions.smoother = ckfOptions.extensions.smoother; + // extensions.outlierFinder = ckfOptions.extensions.outlierFinder; + + KalmanFitterOptions kfOptions( + ckfOptions.geoContext, ckfOptions.magFieldContext, + ckfOptions.calibrationContext, extensions, + ckfOptions.propagatorPlainOptions, ckfOptions.referenceSurface, + ckfOptions.multipleScattering, ckfOptions.energyLoss, false); + + return m_kalmanFitter.fit(it, end, sParameters, kfOptions, trackContainer); + } + + template class holder_t> + auto findTracks( + sli_kf_t it, sli_kf_t end, const start_parameters_t& kfStartParameters, + const CombinatorialKalmanFilterOptions& ckfOptions, + TrackContainer& trackContainer) const + -> Result> { + // The KF run + auto kfResult = + runKalmanFitter(it, end, kfStartParameters, ckfOptions, trackContainer); + + if (!kfResult.ok()) { + return kfResult.error(); + } + + ACTS_INFO("Done KF fitting"); + + // The CKF run + using SourceLinkAccessor = SourceLinkAccessorDelegate; + + using ThisCkf = CombinatorialKalmanFilter; + using Aborter = typename ThisCkf::template Aborter; + using Actor = typename ThisCkf::template Actor; + using Actors = ActionList; + using Aborters = AbortList; + + PropagatorOptions propOptions(ckfOptions.geoContext, + ckfOptions.magFieldContext); + + propOptions.setPlainOptions(ckfOptions.propagatorPlainOptions); + + auto& combKalmanActor = propOptions.actionList.template get(); + combKalmanActor.targetSurface = ckfOptions.referenceSurface; + combKalmanActor.multipleScattering = ckfOptions.multipleScattering; + combKalmanActor.energyLoss = ckfOptions.energyLoss; + combKalmanActor.smoothing = ckfOptions.smoothing; + combKalmanActor.m_sourcelinkAccessor = ckfOptions.sourcelinkAccessor; + combKalmanActor.m_extensions = ckfOptions.extensions; + + // Prepare the start parameters + const auto state = *kfResult->trackStates().end(); + + BoundTrackParameters ckfStartParameters( + state.referenceSurface().getSharedPtr(), state.filtered(), + state.filteredCovariance()); + + // Prepare the result + using CkfResultType = + typename propagator_t::template action_list_t_result_t< + CurvilinearTrackParameters, Actors>; + CkfResultType inputResult; + + auto& ckfResult = + inputResult.template get>(); + ckfResult.fittedStates = &kfResult->container().trackStateContainer(); + ckfResult.lastMeasurementIndices.push_back(kfResult->tipIndex()); + ckfResult.lastTrackIndices.push_back(kfResult->tipIndex()); + + const CombinatorialKalmanFilterTipState tipState{ + kfResult->nMeasurements() + kfResult->nHoles(), + kfResult->nTrackStates(), kfResult->nMeasurements(), 0ul, + kfResult->nHoles()}; + + ckfResult.activeTips.push_back({kfResult->tipIndex(), tipState}); + + // Run the CombinatorialKalmanFilter. + auto result = m_propagator.template propagate( + ckfStartParameters, propOptions, std::move(inputResult)); + + if (!result.ok()) { + ACTS_ERROR("Propapation failed: " << result.error() << " " + << result.error().message() + << " with the initial parameters:\n" + << ckfStartParameters.parameters()); + return result.error(); + } + + ACTS_INFO("Done CKF fitting"); + + return result->template get>(); + } +}; // namespace Acts + +} // namespace Acts diff --git a/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFindingX/ParameterFromTrajectoryAlgorithm.hpp b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFindingX/ParameterFromTrajectoryAlgorithm.hpp new file mode 100644 index 00000000000..f809c58a005 --- /dev/null +++ b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFindingX/ParameterFromTrajectoryAlgorithm.hpp @@ -0,0 +1,73 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2022 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/EventData/MultiTrajectory.hpp" +#include "ActsExamples/EventData/Trajectories.hpp" +#include "ActsExamples/Framework/BareAlgorithm.hpp" +#include "ActsExamples/Framework/WhiteBoard.hpp" + +#include +#include + +namespace ActsExamples { + +class ParameterFromTrajectoryAlgorithm final : public BareAlgorithm { + public: + struct Config { + /// Input spacepoints collection. + std::string inputTrajectories; + + /// Output protoTracks collection. + std::string outputParamters; + }; + + /// Constructor of the track finding algorithm + /// + /// @param cfg is the config struct to configure the algorithm + /// @param level is the logging level + ParameterFromTrajectoryAlgorithm(Config cfg, Acts::Logging::Level lvl) + : BareAlgorithm("ParsFromTraj", lvl), m_cfg(cfg) {} + + virtual ~ParameterFromTrajectoryAlgorithm() {} + + /// Filter the measurements + /// + /// @param ctx is the algorithm context that holds event-wise information + /// @return a process code to steer the algorithm flow + ActsExamples::ProcessCode execute( + const ActsExamples::AlgorithmContext& ctx) const final { + const auto& trajs = + ctx.eventStore.get(m_cfg.inputTrajectories); + + TrackParametersContainer trackParameters; + + for (const auto& traj : trajs) { + const auto i = traj.tips().front(); + const auto state = traj.multiTrajectory().getTrackState(i); + + trackParameters.emplace_back(state.referenceSurface().getSharedPtr(), + state.smoothed(), + state.smoothedCovariance()); + } + + ctx.eventStore.add(m_cfg.outputParamters, + std::move(trackParameters)); + + return ProcessCode::SUCCESS; + } + + const Config& config() const { return m_cfg; } + + private: + // configuration + Config m_cfg; +}; + +} // namespace ActsExamples diff --git a/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFindingX/SourceLinkSelectorAlgorithm.hpp b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFindingX/SourceLinkSelectorAlgorithm.hpp new file mode 100644 index 00000000000..4b1ca8ef0d1 --- /dev/null +++ b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFindingX/SourceLinkSelectorAlgorithm.hpp @@ -0,0 +1,78 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2022 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/Geometry/GeometryIdentifier.hpp" +#include "ActsExamples/EventData/IndexSourceLink.hpp" +#include "ActsExamples/Framework/BareAlgorithm.hpp" +#include "ActsExamples/Framework/WhiteBoard.hpp" + +#include +#include + +namespace ActsExamples { + +class SourceLinkSelectorAlgorithm final : public BareAlgorithm { + public: + struct Config { + /// Input spacepoints collection. + std::string inputSourceLinks; + + /// Output protoTracks collection. + std::string outputSourceLinks; + + /// What spacepoints to keep + std::vector geometrySelection; + }; + + /// Constructor of the track finding algorithm + /// + /// @param cfg is the config struct to configure the algorithm + /// @param level is the logging level + SourceLinkSelectorAlgorithm(Config cfg, Acts::Logging::Level lvl) + : BareAlgorithm("SourceLinkSelection", lvl), m_cfg(cfg) {} + + virtual ~SourceLinkSelectorAlgorithm() {} + + /// Filter the measurements + /// + /// @param ctx is the algorithm context that holds event-wise information + /// @return a process code to steer the algorithm flow + ActsExamples::ProcessCode execute( + const ActsExamples::AlgorithmContext& ctx) const final { + const auto& inputSourceLinks = + ctx.eventStore.get(m_cfg.inputSourceLinks); + + IndexSourceLinkContainer outputSourceLinks; + + for (const auto geoId : m_cfg.geometrySelection) { + auto range = selectLowestNonZeroGeometryObject(inputSourceLinks, geoId); + auto groupedByModule = makeGroupBy(range, detail::GeometryIdGetter()); + + for (auto [moduleGeoId, moduleSourceLinks] : groupedByModule) { + for (auto& sourceLink : moduleSourceLinks) { + outputSourceLinks.insert(sourceLink); + } + } + } + + ctx.eventStore.add(m_cfg.outputSourceLinks, + std::move(outputSourceLinks)); + + return ProcessCode::SUCCESS; + } + + const Config& config() const { return m_cfg; } + + private: + // configuration + Config m_cfg; +}; + +} // namespace ActsExamples diff --git a/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp b/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp index bebf8711e20..9579387d303 100644 --- a/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp +++ b/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp @@ -37,6 +37,8 @@ #include #include +#include +#include #include ActsExamples::TrackFindingAlgorithm::TrackFindingAlgorithm( @@ -139,6 +141,12 @@ ActsExamples::ProcessCode ActsExamples::TrackFindingAlgorithm::execute( unsigned int nSeed = 0; + std::vector nTracksPerSeeds; + nTracksPerSeeds.reserve(initialParameters.size()); + + std::vector nSelTracksPerSeeds; + nSelTracksPerSeeds.reserve(initialParameters.size()); + for (std::size_t iseed = 0; iseed < initialParameters.size(); ++iseed) { // Clear trackContainerTemp and trackStateContainerTemp tracksTemp.clear(); @@ -156,14 +164,21 @@ ActsExamples::ProcessCode ActsExamples::TrackFindingAlgorithm::execute( } auto& tracksForSeed = result.value(); + + nTracksPerSeeds.push_back(tracksForSeed.size()); + std::size_t nSelTracks = 0; + for (auto& track : tracksForSeed) { seedNumber(track) = nSeed; if (!m_trackSelector.has_value() || m_trackSelector->isValidTrack(track)) { + nSelTracks++; auto destProxy = tracks.getTrack(tracks.addTrack()); destProxy.copyFrom(track, true); // make sure we copy track states! } } + + nSelTracksPerSeeds.push_back(nSelTracks); } // Compute shared hits from all the reconstructed tracks @@ -171,6 +186,15 @@ ActsExamples::ProcessCode ActsExamples::TrackFindingAlgorithm::execute( computeSharedHits(sourceLinks, tracks); } + { + std::lock_guard guard(m_mutex); + + std::copy(nTracksPerSeeds.begin(), nTracksPerSeeds.end(), + std::back_inserter(m_nTracksPerSeeds)); + std::copy(nSelTracksPerSeeds.begin(), nSelTracksPerSeeds.end(), + std::back_inserter(m_nSelTracksPerSeeds)); + } + ACTS_DEBUG("Finalized track finding with " << tracks.size() << " track candidates."); @@ -198,6 +222,25 @@ ActsExamples::ProcessCode ActsExamples::TrackFindingAlgorithm::finalize() { ACTS_INFO("- failure ratio: " << static_cast(m_nFailedSeeds) / m_nTotalSeeds); + namespace ba = boost::accumulators; + using Accumulator = ba::accumulator_set< + float, ba::features>; + + Accumulator totalAcc; + std::for_each(m_nTracksPerSeeds.begin(), m_nTracksPerSeeds.end(), + [&](auto v) { totalAcc(static_cast(v)); }); + ACTS_INFO("- total number tracks: " << ba::sum(totalAcc)); + ACTS_INFO("- avg tracks per seed: " << ba::mean(totalAcc) << " +- " + << std::sqrt(ba::variance(totalAcc))); + + Accumulator selAcc; + std::for_each(m_nSelTracksPerSeeds.begin(), m_nSelTracksPerSeeds.end(), + [&](auto v) { selAcc(static_cast(v)); }); + ACTS_INFO("- total number tracks (selected only): " << ba::sum(selAcc)); + ACTS_INFO("- avg tracks per seed (selected only): " + << ba::mean(selAcc) << " +- " << std::sqrt(ba::variance(selAcc))); + + // Memory statistics auto memoryStatistics = m_memoryStatistics.combine([](const auto& a, const auto& b) { Acts::VectorMultiTrajectory::Statistics c; diff --git a/Examples/Algorithms/TrackFinding/src/TrackParamsEstimationAlgorithm.cpp b/Examples/Algorithms/TrackFinding/src/TrackParamsEstimationAlgorithm.cpp index 34ab7cd9317..87bdd818e20 100644 --- a/Examples/Algorithms/TrackFinding/src/TrackParamsEstimationAlgorithm.cpp +++ b/Examples/Algorithms/TrackFinding/src/TrackParamsEstimationAlgorithm.cpp @@ -67,7 +67,7 @@ ActsExamples::TrackParamsEstimationAlgorithm::TrackParamsEstimationAlgorithm( ActsExamples::ProcessCode ActsExamples::TrackParamsEstimationAlgorithm::execute( const ActsExamples::AlgorithmContext& ctx) const { auto const& seeds = m_inputSeeds(ctx); - ACTS_VERBOSE("Read " << seeds.size() << " seeds"); + ACTS_DEBUG("Read " << seeds.size() << " seeds"); TrackParametersContainer trackParameters; trackParameters.reserve(seeds.size()); @@ -142,7 +142,7 @@ ActsExamples::ProcessCode ActsExamples::TrackParamsEstimationAlgorithm::execute( } } - ACTS_VERBOSE("Estimated " << trackParameters.size() << " track parameters"); + ACTS_DEBUG("Estimated " << trackParameters.size() << " track parameters"); m_outputTrackParameters(ctx, std::move(trackParameters)); if (m_outputSeeds.isInitialized()) { diff --git a/Examples/Algorithms/TrackFindingExaTrkX/CMakeLists.txt b/Examples/Algorithms/TrackFindingExaTrkX/CMakeLists.txt index b4a0f280555..d66ba5187cf 100644 --- a/Examples/Algorithms/TrackFindingExaTrkX/CMakeLists.txt +++ b/Examples/Algorithms/TrackFindingExaTrkX/CMakeLists.txt @@ -1,6 +1,7 @@ add_library( ActsExamplesTrackFindingExaTrkX SHARED src/TrackFindingAlgorithmExaTrkX.cpp + src/ProtoTrackEffPurPrinter.cpp src/PrototracksToParameters.cpp src/TrackFindingFromPrototrackAlgorithm.cpp ) diff --git a/Examples/Algorithms/TrackFindingExaTrkX/include/ActsExamples/TrackFindingExaTrkX/ProtoTrackEffPurPrinter.hpp b/Examples/Algorithms/TrackFindingExaTrkX/include/ActsExamples/TrackFindingExaTrkX/ProtoTrackEffPurPrinter.hpp new file mode 100644 index 00000000000..4bde2268582 --- /dev/null +++ b/Examples/Algorithms/TrackFindingExaTrkX/include/ActsExamples/TrackFindingExaTrkX/ProtoTrackEffPurPrinter.hpp @@ -0,0 +1,63 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2023 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/EventData/MultiTrajectory.hpp" +#include "ActsExamples/EventData/ProtoTrack.hpp" +#include "ActsExamples/Framework/DataHandle.hpp" +#include "ActsExamples/Framework/IAlgorithm.hpp" +#include "ActsExamples/Framework/WhiteBoard.hpp" +#include "ActsExamples/Utilities/EventDataTransforms.hpp" + +#include +#include + +#include + +namespace ActsExamples { + +class ProtoTrackEffPurPrinter final : public IAlgorithm { + public: + struct Config { + std::string testProtoTracks; + std::string refProtoTracks; + }; + + ProtoTrackEffPurPrinter(Config cfg, Acts::Logging::Level lvl); + + ActsExamples::ProcessCode execute( + const ActsExamples::AlgorithmContext &context) const override; + + ActsExamples::ProcessCode finalize() override; + + const Config &config() const { return m_cfg; } + + private: + Config m_cfg; + + ReadDataHandle m_testProtoTracks{this, + "InputTestProtoTracks"}; + ReadDataHandle m_refProtoTracks{this, + "InputRefProtoTracks"}; + + using Hist = decltype(boost::histogram::make_histogram( + std::declval>())); + + mutable Hist m_effHistogram; + mutable Hist m_purHistogram; + + using CountHist = decltype(boost::histogram::make_histogram( + std::declval>())); + + mutable CountHist m_countPerTrackHist; + + mutable std::mutex m_histogramMutex; +}; + +} // namespace ActsExamples diff --git a/Examples/Algorithms/TrackFindingExaTrkX/include/ActsExamples/TrackFindingExaTrkX/PrototracksToParsAndSeeds.hpp b/Examples/Algorithms/TrackFindingExaTrkX/include/ActsExamples/TrackFindingExaTrkX/PrototracksToParsAndSeeds.hpp new file mode 100644 index 00000000000..f4dab711dd2 --- /dev/null +++ b/Examples/Algorithms/TrackFindingExaTrkX/include/ActsExamples/TrackFindingExaTrkX/PrototracksToParsAndSeeds.hpp @@ -0,0 +1,84 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2023 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +#include "ActsExamples/EventData/ProtoTrack.hpp" +#include "ActsExamples/EventData/SimSeed.hpp" +#include "ActsExamples/EventData/Track.hpp" +#include "ActsExamples/Framework/DataHandle.hpp" +#include "ActsExamples/Framework/IAlgorithm.hpp" + +namespace ActsExamples { + +class PrototracksToParsAndSeeds final : public IAlgorithm { + public: + struct Config { + std::string inputProtoTracks; + std::string inputSpacePoints; + std::string outputSeeds = "seeds-from-prototracks"; + std::string outputProtoTracks = "remaining-prototracks"; + std::string outputParameters = "parameters"; + + // The tracking geometry + std::shared_ptr geometry; + + // Wether to make tight seeds (closest hits to beampipe) or large seeds + bool buildTightSeeds = false; + + /// The minimum magnetic field to trigger the track parameters estimation + double bFieldMin = 0.1 * Acts::UnitConstants::T; + /// Constant term of the loc0 resolution. + double sigmaLoc0 = 25 * Acts::UnitConstants::um; + /// Constant term of the loc1 resolution. + double sigmaLoc1 = 100 * Acts::UnitConstants::um; + /// Phi angular resolution. + double sigmaPhi = 0.02 * Acts::UnitConstants::degree; + /// Theta angular resolution. + double sigmaTheta = 0.02 * Acts::UnitConstants::degree; + /// q/p resolution. + double sigmaQOverP = 0.1 / Acts::UnitConstants::GeV; + /// Time resolution. + double sigmaT0 = 10 * Acts::UnitConstants::ns; + /// Inflate initial covariance. + std::array initialVarInflation = {1., 1., 1., 1., 1., 1.}; + }; + + /// Construct the algorithm. + /// + /// @param cfg is the algorithm configuration + /// @param lvl is the logging level + PrototracksToParsAndSeeds(Config cfg, Acts::Logging::Level lvl); + + ~PrototracksToParsAndSeeds(); + + /// Run the algorithm. + /// + /// @param ctx is the algorithm context with event information + /// @return a process code indication success or failure + ProcessCode execute(const AlgorithmContext& ctx) const final; + + /// Const access to the config + const Config& config() const { return m_cfg; } + + private: + Config m_cfg; + Acts::BoundSquareMatrix m_covariance = Acts::BoundSquareMatrix::Zero(); + + WriteDataHandle m_outputSeeds{this, "OutputSeeds"}; + WriteDataHandle m_outputProtoTracks{this, + "OutputProtoTracks"}; + WriteDataHandle m_outputParameters{ + this, "OutputParameters"}; + ReadDataHandle m_inputSpacePoints{this, + "InputSpacePoints"}; + ReadDataHandle m_inputProtoTracks{this, + "InputProtoTracks"}; +}; + +} // namespace ActsExamples diff --git a/Examples/Algorithms/TrackFindingExaTrkX/include/ActsExamples/TrackFindingExaTrkX/TrackFindingAlgorithmExaTrkX.hpp b/Examples/Algorithms/TrackFindingExaTrkX/include/ActsExamples/TrackFindingExaTrkX/TrackFindingAlgorithmExaTrkX.hpp index d6366d25b43..6e01c7c35dd 100644 --- a/Examples/Algorithms/TrackFindingExaTrkX/include/ActsExamples/TrackFindingExaTrkX/TrackFindingAlgorithmExaTrkX.hpp +++ b/Examples/Algorithms/TrackFindingExaTrkX/include/ActsExamples/TrackFindingExaTrkX/TrackFindingAlgorithmExaTrkX.hpp @@ -11,6 +11,7 @@ #include "Acts/Definitions/Units.hpp" #include "Acts/Plugins/ExaTrkX/ExaTrkXPipeline.hpp" #include "Acts/Plugins/ExaTrkX/Stages.hpp" +#include "Acts/Plugins/ExaTrkX/TorchGraphStoreHook.hpp" #include "ActsExamples/EventData/Cluster.hpp" #include "ActsExamples/EventData/ProtoTrack.hpp" #include "ActsExamples/EventData/SimHit.hpp" @@ -52,6 +53,9 @@ class TrackFindingAlgorithmExaTrkX final : public IAlgorithm { /// Output protoTracks collection. std::string outputProtoTracks; + /// Output graph (optional) + std::string outputGraph; + std::shared_ptr graphConstructor; std::vector> edgeClassifiers; @@ -114,6 +118,8 @@ class TrackFindingAlgorithmExaTrkX final : public IAlgorithm { WriteDataHandle m_outputProtoTracks{this, "OutputProtoTracks"}; + WriteDataHandle m_outputGraph{ + this, "OutputGraph"}; // for truth graph ReadDataHandle m_inputSimHits{this, "InputSimHits"}; diff --git a/Examples/Algorithms/TrackFindingExaTrkX/src/ProtoTrackEffPurPrinter.cpp b/Examples/Algorithms/TrackFindingExaTrkX/src/ProtoTrackEffPurPrinter.cpp new file mode 100644 index 00000000000..fb6e1b55c8c --- /dev/null +++ b/Examples/Algorithms/TrackFindingExaTrkX/src/ProtoTrackEffPurPrinter.cpp @@ -0,0 +1,134 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2023 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "ActsExamples/TrackFindingExaTrkX/ProtoTrackEffPurPrinter.hpp" + +#include +#include + +#include + +namespace bh = boost::histogram; + +ActsExamples::ProtoTrackEffPurPrinter::ProtoTrackEffPurPrinter( + Config cfg, Acts::Logging::Level lvl) + : IAlgorithm("ProtoTrackEfficencyPrinter", lvl), + m_cfg(cfg), + m_effHistogram(bh::make_histogram(bh::axis::regular<>(10, 0.0, 1.0))), + m_purHistogram(bh::make_histogram(bh::axis::regular<>(10, 0.0, 1.0))), + m_countPerTrackHist(bh::make_histogram(bh::axis::integer<>(1, 5))) { + m_testProtoTracks.initialize(m_cfg.testProtoTracks); + m_refProtoTracks.initialize(m_cfg.refProtoTracks); +} + +ActsExamples::ProcessCode ActsExamples::ProtoTrackEffPurPrinter::execute( + const ActsExamples::AlgorithmContext &context) const { + auto testTracks = m_testProtoTracks(context); + auto truthTracks = m_refProtoTracks(context); + + ACTS_INFO("Receiving " << truthTracks.size() << " reference tracks"); + truthTracks.erase(std::remove_if(truthTracks.begin(), truthTracks.end(), + [](const auto &t) { return t.size() < 3; }), + truthTracks.end()); + ACTS_INFO(" -> " << truthTracks.size() << " tracks with size >= 3"); + ACTS_INFO("Receiving " << testTracks.size() << " test tracks"); + + // Build id-to-truth-track map + // For now we assume that each space point only belongs to one truth track + // (this might be wrong if we merge clusters) + std::unordered_map idToTruthTrack; + + for (auto i = 0ul; i < truthTracks.size(); ++i) { + for (auto el : truthTracks[i]) { + idToTruthTrack[el] = i; + } + } + + // Go through test tracks and search truth track ids + constexpr static std::size_t invalid = + std::numeric_limits::max(); + + std::vector testTrackPurities(testTracks.size(), 0.f); + std::vector trueTrackEfficiencies(truthTracks.size(), 0.f); + std::vector trueTracksPerTestTrack; + + { + // allocate once for memory optimization + std::vector truthTrackIds; + std::vector truthTrackUniqueIds; + + for (auto testId = 0ul; testId < testTracks.size(); ++testId) { + const auto &testTrack = testTracks[testId]; + + // Build vector of truth track ids + for (const auto el : testTrack) { + // invalid means, the points is not associated to any track (noise) + auto tid = idToTruthTrack.count(el) > 0 ? idToTruthTrack[el] : invalid; + truthTrackIds.push_back(tid); + } + + // Find out which truth track ids are most fequent + std::sort(truthTrackIds.begin(), truthTrackIds.end()); + std::unique_copy(truthTrackIds.begin(), truthTrackIds.end(), + std::back_inserter(truthTrackUniqueIds)); + + std::sort( + truthTrackUniqueIds.begin(), truthTrackUniqueIds.end(), + [&](const auto &a, const auto &b) { + auto ac = std::count(truthTrackIds.begin(), truthTrackIds.end(), a); + auto bc = std::count(truthTrackIds.begin(), truthTrackIds.end(), b); + return ac > bc; // sort many-to-few + }); + + trueTracksPerTestTrack.push_back(truthTrackUniqueIds.size()); + + // compute metrics + if (truthTrackUniqueIds[0] == invalid && + truthTrackUniqueIds.size() == 1) { + testTrackPurities[testId] = 0.f; + } else { + const auto truthId = truthTrackUniqueIds[0] == invalid + ? truthTrackUniqueIds[1] + : truthTrackUniqueIds[0]; + const auto nhits = + std::count(truthTrackIds.begin(), truthTrackIds.end(), truthId); + + const auto &truthTrack = truthTracks[truthId]; + + // Ensure eff,pur < 1.0, so the binning in the histogram is nice + float eff = static_cast(nhits) / truthTrack.size(); + eff = std::min(eff, 0.9999f); + + float pur = static_cast(nhits) / testTrack.size(); + pur = std::min(pur, 0.9999f); + + trueTrackEfficiencies[truthId] = + std::max(trueTrackEfficiencies[truthId], eff); + testTrackPurities[testId] = pur; + } + + // clear vectors + truthTrackIds.clear(); + truthTrackUniqueIds.clear(); + } + } + + std::lock_guard{m_histogramMutex}; + m_effHistogram.fill(trueTrackEfficiencies); + m_purHistogram.fill(testTrackPurities); + m_countPerTrackHist.fill(trueTracksPerTestTrack); + + return {}; +} + +ActsExamples::ProcessCode ActsExamples::ProtoTrackEffPurPrinter::finalize() { + ACTS_INFO("Truth track efficiency:\n" << m_effHistogram); + ACTS_INFO("Test track purity:\n" << m_purHistogram); + ACTS_INFO("Particles per test track:\n" << m_countPerTrackHist); + return {}; +} diff --git a/Examples/Algorithms/TrackFindingExaTrkX/src/PrototracksToParsAndSeeds.cpp b/Examples/Algorithms/TrackFindingExaTrkX/src/PrototracksToParsAndSeeds.cpp new file mode 100644 index 00000000000..75dc4dfc4db --- /dev/null +++ b/Examples/Algorithms/TrackFindingExaTrkX/src/PrototracksToParsAndSeeds.cpp @@ -0,0 +1,203 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2023 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "ActsExamples/TrackFindingExaTrkX/PrototracksToParsAndSeeds.hpp" + +#include "Acts/Seeding/BinFinder.hpp" +#include "Acts/Seeding/BinnedSPGroup.hpp" +#include "Acts/Seeding/EstimateTrackParamsFromSeed.hpp" +#include "Acts/Seeding/InternalSpacePoint.hpp" +#include "Acts/Seeding/SeedFilter.hpp" +#include "Acts/Seeding/SeedFinder.hpp" +#include "Acts/Seeding/SeedFinderConfig.hpp" +#include "Acts/Utilities/Zip.hpp" +#include "ActsExamples/EventData/IndexSourceLink.hpp" +#include "ActsExamples/EventData/ProtoTrack.hpp" +#include "ActsExamples/EventData/SimSeed.hpp" +#include "ActsExamples/Framework/WhiteBoard.hpp" +#include "ActsExamples/Utilities/EventDataTransforms.hpp" + +#include + +using namespace ActsExamples; +using namespace Acts::UnitLiterals; + +namespace ActsExamples { + +PrototracksToParsAndSeeds::PrototracksToParsAndSeeds(Config cfg, + Acts::Logging::Level lvl) + : IAlgorithm("PrototracksToParsAndSeeds", lvl), m_cfg(std::move(cfg)) { + m_outputSeeds.initialize(m_cfg.outputSeeds); + m_outputProtoTracks.initialize(m_cfg.outputProtoTracks); + m_inputProtoTracks.initialize(m_cfg.inputProtoTracks); + m_inputSpacePoints.initialize(m_cfg.inputSpacePoints); + m_outputParameters.initialize(m_cfg.outputParameters); + + if (m_cfg.geometry == nullptr) { + throw std::invalid_argument("No geometry given"); + } + + // m_advancedSeeding = std::make_unique(logger()); + + // Set up the track parameters covariance (the same for all tracks) + m_covariance(Acts::eBoundLoc0, Acts::eBoundLoc0) = + m_cfg.initialVarInflation[Acts::eBoundLoc0] * cfg.sigmaLoc0 * + m_cfg.sigmaLoc0; + m_covariance(Acts::eBoundLoc1, Acts::eBoundLoc1) = + m_cfg.initialVarInflation[Acts::eBoundLoc1] * cfg.sigmaLoc1 * + m_cfg.sigmaLoc1; + m_covariance(Acts::eBoundPhi, Acts::eBoundPhi) = + m_cfg.initialVarInflation[Acts::eBoundPhi] * cfg.sigmaPhi * + m_cfg.sigmaPhi; + m_covariance(Acts::eBoundTheta, Acts::eBoundTheta) = + m_cfg.initialVarInflation[Acts::eBoundTheta] * cfg.sigmaTheta * + m_cfg.sigmaTheta; + m_covariance(Acts::eBoundQOverP, Acts::eBoundQOverP) = + m_cfg.initialVarInflation[Acts::eBoundQOverP] * cfg.sigmaQOverP * + m_cfg.sigmaQOverP; + m_covariance(Acts::eBoundTime, Acts::eBoundTime) = + m_cfg.initialVarInflation[Acts::eBoundTime] * m_cfg.sigmaT0 * + m_cfg.sigmaT0; +} + +PrototracksToParsAndSeeds::~PrototracksToParsAndSeeds() {} + +ProcessCode PrototracksToParsAndSeeds::execute( + const AlgorithmContext &ctx) const { + const auto &sps = m_inputSpacePoints(ctx); + auto prototracks = m_inputProtoTracks(ctx); + + // Make some lookup tables. Allocate space for the maximum number of indices + // (max 2 source links per spacepoint) + std::vector indexToSpacepoint(2 * sps.size(), nullptr); + std::vector indexToGeoId( + 2 * sps.size(), Acts::GeometryIdentifier{0}); + + for (const auto &sp : sps) { + for (const auto &sl : sp.sourceLinks()) { + const auto &isl = sl.template get(); + indexToSpacepoint[isl.index()] = &sp; + indexToGeoId[isl.index()] = isl.geometryId(); + } + } + + ProtoTrackContainer seededTracks; + seededTracks.reserve(prototracks.size()); + + SimSeedContainer seeds; + seeds.reserve(prototracks.size()); + + TrackParametersContainer parameters; + parameters.reserve(prototracks.size()); + + // Loop over the prototracks to make seeds + ProtoTrack tmpTrack; + std::vector tmpSps; + std::size_t skippedTracks = 0; + for (auto &track : prototracks) { + ACTS_VERBOSE("Try to get seed from prototrack with " << track.size() + << " hits"); + // Make prototrack unique with respect to volume and layer + // so we don't get a seed where we have two spacepoints on the same layer + + // Here, we want to create a seed only if the prototrack with removed unique + // layer-volume spacepoints has 3 or more hits. However, if this is the + // case, we want to keep the whole prototrack. Therefore, we operate on a + // tmpTrack. + std::sort(track.begin(), track.end(), [&](auto a, auto b) { + if (indexToGeoId[a].volume() != indexToGeoId[b].volume()) { + return indexToGeoId[a].volume() < indexToGeoId[b].volume(); + } + return indexToGeoId[a].layer() < indexToGeoId[b].layer(); + }); + + tmpTrack.clear(); + std::unique_copy( + track.begin(), track.end(), std::back_inserter(tmpTrack), + [&](auto a, auto b) { + return indexToGeoId[a].volume() == indexToGeoId[b].volume() && + indexToGeoId[a].layer() == indexToGeoId[b].layer(); + }); + + // in this case we cannot seed properly + if (tmpTrack.size() < 3) { + ACTS_DEBUG( + "Cannot seed because less then three hits with unique (layer, " + "volume)"); + skippedTracks++; + continue; + } + + // Make the seed + tmpSps.clear(); + std::transform(track.begin(), track.end(), std::back_inserter(tmpSps), + [&](auto i) { return indexToSpacepoint[i]; }); + tmpSps.erase(std::remove_if(tmpSps.begin(), tmpSps.end(), + [](auto sp) { return sp == nullptr; }), + tmpSps.end()); + + if (tmpSps.size() < 3) { + ACTS_WARNING("Could not find all spacepoints, skip"); + skippedTracks++; + continue; + } + + std::sort(tmpSps.begin(), tmpSps.end(), + [](const auto &a, const auto &b) { return a->r() < b->r(); }); + + // Simply use r = m*z + t and solve for r=0 to find z vertex position... + // Probably not the textbook way to do + const auto m = (tmpSps.back()->r() - tmpSps.front()->r()) / + (tmpSps.back()->z() - tmpSps.front()->z()); + const auto t = tmpSps.front()->r() - m * tmpSps.front()->z(); + const auto z_vertex = -t / m; + const auto s = tmpSps.size(); + + SimSeed seed = + m_cfg.buildTightSeeds + ? SimSeed(*tmpSps[0], *tmpSps[1], *tmpSps[2], z_vertex) + : SimSeed(*tmpSps[0], *tmpSps[s / 2], *tmpSps[s - 1], z_vertex); + + // Compute parameters + const auto geoId = seed.sp() + .front() + ->sourceLinks() + .front() + .template get() + .geometryId(); + const auto &surface = *m_cfg.geometry->findSurface(geoId); + + auto pars = Acts::estimateTrackParamsFromSeed( + {}, seed.sp().begin(), seed.sp().end(), surface, {0., 0., 2_T}, 0.0); + + if (not pars) { + ACTS_WARNING("Skip track because of bad params"); + } + + seededTracks.push_back(track); + seeds.emplace_back(std::move(seed)); + parameters.push_back( + Acts::BoundTrackParameters(surface.getSharedPtr(), *pars, m_covariance, + Acts::ParticleHypothesis::pion())); + } + + if (skippedTracks > 0) { + ACTS_WARNING("Skipped seeding of " << skippedTracks); + } + + ACTS_DEBUG("Seeded " << seeds.size() << " out of " << prototracks.size() + << " prototracks"); + + m_outputSeeds(ctx, std::move(seeds)); + m_outputProtoTracks(ctx, std::move(seededTracks)); + m_outputParameters(ctx, std::move(parameters)); + + return ProcessCode::SUCCESS; +} + +} // namespace ActsExamples diff --git a/Examples/Algorithms/TrackFindingExaTrkX/src/TrackFindingAlgorithmExaTrkX.cpp b/Examples/Algorithms/TrackFindingExaTrkX/src/TrackFindingAlgorithmExaTrkX.cpp index 4215828d62c..4bf1dc91b7b 100644 --- a/Examples/Algorithms/TrackFindingExaTrkX/src/TrackFindingAlgorithmExaTrkX.cpp +++ b/Examples/Algorithms/TrackFindingExaTrkX/src/TrackFindingAlgorithmExaTrkX.cpp @@ -9,6 +9,7 @@ #include "ActsExamples/TrackFindingExaTrkX/TrackFindingAlgorithmExaTrkX.hpp" #include "Acts/Definitions/Units.hpp" +#include "Acts/Plugins/ExaTrkX/TorchGraphStoreHook.hpp" #include "Acts/Plugins/ExaTrkX/TorchTruthGraphMetricsHook.hpp" #include "Acts/Utilities/Zip.hpp" #include "ActsExamples/EventData/Index.hpp" @@ -17,6 +18,7 @@ #include "ActsExamples/EventData/SimSpacePoint.hpp" #include "ActsExamples/Framework/WhiteBoard.hpp" +#include #include using namespace ActsExamples; @@ -31,6 +33,7 @@ class ExamplesEdmHook : public Acts::ExaTrkXHook { std::unique_ptr m_logger; std::unique_ptr m_truthGraphHook; std::unique_ptr m_targetGraphHook; + std::unique_ptr m_graphStoreHook; const Acts::Logger& logger() const { return *m_logger; } @@ -70,6 +73,8 @@ class ExamplesEdmHook : public Acts::ExaTrkXHook { std::vector truthGraph; std::vector targetGraph; + std::size_t notMatched = 0; + for (auto& [pid, track] : tracks) { // Sort by hit index, so the edges are connected correctly std::sort(track.begin(), track.end(), [](const auto& a, const auto& b) { @@ -78,7 +83,8 @@ class ExamplesEdmHook : public Acts::ExaTrkXHook { auto found = particles.find(pid); if (found == particles.end()) { - ACTS_WARNING("Did not find " << pid << ", skip track"); + ACTS_VERBOSE("Did not find " << pid << ", cannot add to target graph"); + notMatched++; continue; } @@ -94,21 +100,30 @@ class ExamplesEdmHook : public Acts::ExaTrkXHook { } } + ACTS_DEBUG("Was not able to match " + << notMatched + << " particles, these might be missing in target graph"); + m_truthGraphHook = std::make_unique( truthGraph, logger.clone()); m_targetGraphHook = std::make_unique( targetGraph, logger.clone()); + m_graphStoreHook = std::make_unique(); } ~ExamplesEdmHook() {} - void operator()(const std::any& nodes, const std::any& edges) const override { + auto storedGraph() const { return m_graphStoreHook->storedGraph(); } + + void operator()(const std::any& nodes, const std::any& edges, + const std::any& weights) const override { ACTS_INFO("Metrics for total graph:"); - (*m_truthGraphHook)(nodes, edges); + (*m_truthGraphHook)(nodes, edges, weights); ACTS_INFO("Metrics for target graph (pT > " << m_targetPT / Acts::UnitConstants::GeV << " GeV, nHits >= " << m_targetSize << "):"); - (*m_targetGraphHook)(nodes, edges); + (*m_targetGraphHook)(nodes, edges, weights); + (*m_graphStoreHook)(nodes, edges, weights); } }; @@ -153,6 +168,8 @@ ActsExamples::TrackFindingAlgorithmExaTrkX::TrackFindingAlgorithmExaTrkX( m_inputParticles.maybeInitialize(m_cfg.inputParticles); m_inputMeasurementMap.maybeInitialize(m_cfg.inputMeasurementSimhitsMap); + m_outputGraph.maybeInitialize(m_cfg.outputGraph); + // reserve space for timing m_timing.classifierTimes.resize( m_cfg.edgeClassifiers.size(), @@ -174,6 +191,12 @@ ActsExamples::ProcessCode ActsExamples::TrackFindingAlgorithmExaTrkX::execute( const ActsExamples::AlgorithmContext& ctx) const { // Read input data auto spacepoints = m_inputSpacePoints(ctx); +#if 0 + std::sort(spacepoints.begin(), spacepoints.end(), + [](const auto& a, const auto& b) { + return std::hypot(a.x(), a.y()) < std::hypot(b.x(), b.y()); + }); +#endif auto hook = std::make_unique(); if (m_inputSimHits.isInitialized() && m_inputMeasurementMap.isInitialized()) { @@ -239,6 +262,22 @@ ActsExamples::ProcessCode ActsExamples::TrackFindingAlgorithmExaTrkX::execute( ACTS_DEBUG("Avg cell count: " << sumCells / spacepoints.size()); ACTS_DEBUG("Avg activation: " << sumActivation / sumCells); +#if 0 + { + std::stringstream ss; + std::copy(features.begin(), features.begin() + numFeatures, + std::ostream_iterator(ss, " ")); + ss << "\n"; + std::copy(features.end() - numFeatures, features.end(), + std::ostream_iterator(ss, " ")); + ss << "\n"; + ACTS_DEBUG("First & last row:\n" << ss.str()); + } +#endif + + ACTS_DEBUG("Avg cell count: " << sumCells / spacepoints.size()); + ACTS_DEBUG("Avg activation: " << sumActivation / sumCells); + // Run the pipeline const auto trackCandidates = [&]() { const int deviceHint = -1; @@ -267,15 +306,35 @@ ActsExamples::ProcessCode ActsExamples::TrackFindingAlgorithmExaTrkX::execute( // Make the prototracks std::vector protoTracks; protoTracks.reserve(trackCandidates.size()); + + int nShortTracks = 0; + for (auto& x : trackCandidates) { + if (x.size() < 3) { + nShortTracks++; + continue; + } + ProtoTrack onetrack; + onetrack.reserve(x.size()); + std::copy(x.begin(), x.end(), std::back_inserter(onetrack)); protoTracks.push_back(std::move(onetrack)); } + ACTS_INFO("Removed " << nShortTracks << " with less then 3 hits"); ACTS_INFO("Created " << protoTracks.size() << " proto tracks"); m_outputProtoTracks(ctx, std::move(protoTracks)); + if (auto dhook = dynamic_cast(&*hook); + dhook && m_outputGraph.isInitialized()) { + auto graph = dhook->storedGraph(); + std::transform( + graph.first.begin(), graph.first.end(), graph.first.begin(), + [&](const auto& a) -> int64_t { return spacepointIDs.at(a); }); + m_outputGraph(ctx, std::move(graph)); + } + return ActsExamples::ProcessCode::SUCCESS; } diff --git a/Examples/Algorithms/TrackFindingExaTrkX/src/TrackFindingFromPrototrackAlgorithm.cpp b/Examples/Algorithms/TrackFindingExaTrkX/src/TrackFindingFromPrototrackAlgorithm.cpp index d62ab4de413..8dcb2672fdd 100644 --- a/Examples/Algorithms/TrackFindingExaTrkX/src/TrackFindingFromPrototrackAlgorithm.cpp +++ b/Examples/Algorithms/TrackFindingExaTrkX/src/TrackFindingFromPrototrackAlgorithm.cpp @@ -178,8 +178,6 @@ ActsExamples::ProcessCode TrackFindingFromPrototrackAlgorithm::execute( // once this is done. // Compute shared hits from all the reconstructed tracks if // (m_cfg.computeSharedHits) { - // computeSharedHits(sourceLinks, tracks); - // } ACTS_INFO("Event " << ctx.eventNumber << ": " << nFailed << " / " << nSeed << " failed (" << ((100.f * nFailed) / nSeed) << "%)"); diff --git a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSelector.cpp b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSelector.cpp index 5858f41fc57..74ac8ae2f6f 100644 --- a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSelector.cpp +++ b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSelector.cpp @@ -54,10 +54,10 @@ ActsExamples::ParticleSelector::ParticleSelector(const Config& config, ACTS_DEBUG("remove secondary particles " << m_cfg.removeSecondaries); // We only initialize this if we actually select on this - if (m_cfg.measurementsMin > 0 || + if (m_cfg.measurementsMin > 0 or m_cfg.measurementsMax < std::numeric_limits::max()) { m_inputMap.initialize(m_cfg.inputMeasurementParticlesMap); - ACTS_DEBUG("selection particle number of measurements [" + ACTS_DEBUG("selection particle number of measurments [" << m_cfg.measurementsMin << "," << m_cfg.measurementsMax << ")"); } } @@ -81,7 +81,7 @@ ActsExamples::ProcessCode ActsExamples::ParticleSelector::execute( // helper functions to select tracks auto within = [](auto x, auto min, auto max) { - return (min <= x) && (x < max); + return (min <= x) and (x < max); }; auto isValidParticle = [&](const ActsFatras::Particle& p) { @@ -89,14 +89,14 @@ ActsExamples::ProcessCode ActsExamples::ParticleSelector::execute( const auto phi = Acts::VectorHelpers::phi(p.direction()); const auto rho = Acts::VectorHelpers::perp(p.position()); // define charge selection - const bool validNeutral = (p.charge() == 0) && !m_cfg.removeNeutral; - const bool validCharged = (p.charge() != 0) && !m_cfg.removeCharged; - const bool validCharge = validNeutral || validCharged; - const bool validSecondary = !m_cfg.removeSecondaries || !p.isSecondary(); + const bool validNeutral = (p.charge() == 0) and not m_cfg.removeNeutral; + const bool validCharged = (p.charge() != 0) and not m_cfg.removeCharged; + const bool validCharge = validNeutral or validCharged; + const bool validSecondary = not m_cfg.removeSecondaries or !p.isSecondary(); - nInvalidCharge += static_cast(!validCharge); + nInvalidCharge += not validCharge; - // default valid measurement count to true and only change if we have loaded + // default valid measurment count to true and only change if we have loaded // the measurement particles map bool validMeasurementCount = true; if (particlesMeasMap) { @@ -109,18 +109,17 @@ ActsExamples::ProcessCode ActsExamples::ParticleSelector::execute( << p.particleId()); } - nInvalidMeasurementCount += - static_cast(!validMeasurementCount); + nInvalidMeasurementCount += not validMeasurementCount; - return validCharge && validSecondary && validMeasurementCount && - within(p.transverseMomentum(), m_cfg.ptMin, m_cfg.ptMax) && - within(std::abs(eta), m_cfg.absEtaMin, m_cfg.absEtaMax) && - within(eta, m_cfg.etaMin, m_cfg.etaMax) && - within(phi, m_cfg.phiMin, m_cfg.phiMax) && + return validCharge and validSecondary and validMeasurementCount and + within(p.transverseMomentum(), m_cfg.ptMin, m_cfg.ptMax) and + within(std::abs(eta), m_cfg.absEtaMin, m_cfg.absEtaMax) and + within(eta, m_cfg.etaMin, m_cfg.etaMax) and + within(phi, m_cfg.phiMin, m_cfg.phiMax) and within(std::abs(p.position()[Acts::ePos2]), m_cfg.absZMin, - m_cfg.absZMax) && - within(rho, m_cfg.rhoMin, m_cfg.rhoMax) && - within(p.time(), m_cfg.timeMin, m_cfg.timeMax) && + m_cfg.absZMax) and + within(rho, m_cfg.rhoMin, m_cfg.rhoMax) and + within(p.time(), m_cfg.timeMin, m_cfg.timeMax) and within(p.mass(), m_cfg.mMin, m_cfg.mMax); }; diff --git a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSelector.hpp b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSelector.hpp index d6c9fba3bc9..a7ecab1100b 100644 --- a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSelector.hpp +++ b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSelector.hpp @@ -35,6 +35,8 @@ class ParticleSelector final : public IAlgorithm { std::string inputParticles; /// Input measurement particles map (Optional) std::string inputMeasurementParticlesMap; + /// Input measurements (Optional) + std::string inputMeasurements; /// The output particles collection. std::string outputParticles; // Minimum/maximum distance from the origin in the transverse plane. @@ -84,6 +86,8 @@ class ParticleSelector final : public IAlgorithm { ReadDataHandle m_inputParticles{this, "InputParticles"}; ReadDataHandle> m_inputMap{ this, "InputMeasurementParticlesMap"}; + ReadDataHandle m_inputMeasurements{this, + "InputMeasurements"}; WriteDataHandle m_outputParticles{this, "OutputParticles"}; diff --git a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthTrackFinder.cpp b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthTrackFinder.cpp index f738ba4c49f..9e1afb2e07e 100644 --- a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthTrackFinder.cpp +++ b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthTrackFinder.cpp @@ -62,6 +62,10 @@ ProcessCode TruthTrackFinder::execute(const AlgorithmContext& ctx) const { const auto& hits = makeRange(particleHitsMap.equal_range(particle.particleId())); ACTS_VERBOSE(" - Prototrack from " << hits.size() << " hits"); + if (hits.size() < m_cfg.minHits) { + ACTS_VERBOSE(" --> skip"); + continue; + } // fill hit indices to create the proto track ProtoTrack track; track.reserve(hits.size()); diff --git a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthTrackFinder.hpp b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthTrackFinder.hpp index 8c68906dcf1..51b7f6f04c5 100644 --- a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthTrackFinder.hpp +++ b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/TruthTrackFinder.hpp @@ -43,6 +43,8 @@ class TruthTrackFinder final : public IAlgorithm { std::string inputMeasurementParticlesMap; /// The output proto tracks collection. std::string outputProtoTracks; + /// Minimum Hits + std::size_t minHits = 0; }; TruthTrackFinder(const Config& config, Acts::Logging::Level level); diff --git a/Examples/Algorithms/Utilities/CMakeLists.txt b/Examples/Algorithms/Utilities/CMakeLists.txt index c99737e1305..c9d89d34d18 100644 --- a/Examples/Algorithms/Utilities/CMakeLists.txt +++ b/Examples/Algorithms/Utilities/CMakeLists.txt @@ -5,6 +5,8 @@ add_library( src/TrajectoriesToPrototracks.cpp src/TrackSelectorAlgorithm.cpp src/TracksToTrajectories.cpp + src/HitSelector.cpp + src/MakeMeasurementParticlesMap.cpp src/TracksToParameters.cpp) target_include_directories( ActsExamplesUtilities diff --git a/Examples/Algorithms/Utilities/include/ActsExamples/Utilities/HitSelector.hpp b/Examples/Algorithms/Utilities/include/ActsExamples/Utilities/HitSelector.hpp new file mode 100644 index 00000000000..ee356445670 --- /dev/null +++ b/Examples/Algorithms/Utilities/include/ActsExamples/Utilities/HitSelector.hpp @@ -0,0 +1,50 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2019-2023 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/TrackFinding/TrackSelector.hpp" +#include "Acts/Utilities/Logger.hpp" +#include "ActsExamples/EventData/SimHit.hpp" +#include "ActsExamples/Framework/DataHandle.hpp" +#include "ActsExamples/Framework/IAlgorithm.hpp" + +#include +#include +#include + +namespace ActsExamples { + +/// Select tracks by applying some selection cuts. +class HitSelector final : public IAlgorithm { + public: + struct Config { + /// Input track collection. + std::string inputHits; + /// Output track collection + std::string outputHits; + + /// Time cut + double maxTime = std::numeric_limits::max(); + }; + + HitSelector(const Config& config, Acts::Logging::Level level); + + ProcessCode execute(const AlgorithmContext& ctx) const final; + + /// Get readonly access to the config parameters + const Config& config() const { return m_cfg; } + + private: + Config m_cfg; + + ReadDataHandle m_inputHits{this, "InputHits"}; + WriteDataHandle m_outputHits{this, "OutputHits"}; +}; + +} // namespace ActsExamples diff --git a/Examples/Algorithms/Utilities/include/ActsExamples/Utilities/MakeMeasurementParticlesMap.hpp b/Examples/Algorithms/Utilities/include/ActsExamples/Utilities/MakeMeasurementParticlesMap.hpp new file mode 100644 index 00000000000..eaacaf4bd79 --- /dev/null +++ b/Examples/Algorithms/Utilities/include/ActsExamples/Utilities/MakeMeasurementParticlesMap.hpp @@ -0,0 +1,52 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2023 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +#include "ActsExamples/EventData/Index.hpp" +#include "ActsExamples/EventData/Measurement.hpp" +#include "ActsExamples/EventData/SimHit.hpp" +#include "ActsExamples/Framework/DataHandle.hpp" +#include "ActsExamples/Framework/IAlgorithm.hpp" + +namespace ActsExamples { + +class MakeMeasurementParticlesMap final : public IAlgorithm { + public: + struct Config { + std::string inputSimHits; + std::string inputMeasurementSimhitMap; + std::string outputMeasurementParticlesMap; + }; + + /// Construct the algorithm. + /// + /// @param cfg is the algorithm configuration + /// @param lvl is the logging level + MakeMeasurementParticlesMap(Config cfg, Acts::Logging::Level lvl); + + /// Run the algorithm. + /// + /// @param ctx is the algorithm context with event information + /// @return a process code indication success or failure + ProcessCode execute(const AlgorithmContext& ctx) const final; + + /// Const access to the config + const Config& config() const { return m_cfg; } + + private: + Config m_cfg; + + WriteDataHandle> m_outputParticleMap{ + this, "OutputMeasurementParticlesMap"}; + ReadDataHandle> m_inputHitMap{ + this, "InputMeasurementSimhitMap"}; + ReadDataHandle m_inputHits{this, "InputHits"}; +}; + +} // namespace ActsExamples diff --git a/Examples/Algorithms/Utilities/src/HitSelector.cpp b/Examples/Algorithms/Utilities/src/HitSelector.cpp new file mode 100644 index 00000000000..1da7e059483 --- /dev/null +++ b/Examples/Algorithms/Utilities/src/HitSelector.cpp @@ -0,0 +1,33 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2019-2023 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "ActsExamples/Utilities/HitSelector.hpp" + +ActsExamples::HitSelector::HitSelector(const Config& config, + Acts::Logging::Level level) + : IAlgorithm("HitSelector", level), m_cfg(config) { + m_inputHits.initialize(m_cfg.inputHits); + m_outputHits.initialize(m_cfg.outputHits); +} + +ActsExamples::ProcessCode ActsExamples::HitSelector::execute( + const ActsExamples::AlgorithmContext& ctx) const { + const auto& hits = m_inputHits(ctx); + SimHitContainer selectedHits; + + std::copy_if(hits.begin(), hits.end(), + std::inserter(selectedHits, selectedHits.begin()), + [&](const auto& hit) { return hit.time() < m_cfg.maxTime; }); + + ACTS_DEBUG("selected " << selectedHits.size() << " from " << hits.size() + << " hits"); + + m_outputHits(ctx, std::move(selectedHits)); + + return {}; +} diff --git a/Examples/Algorithms/Utilities/src/MakeMeasurementParticlesMap.cpp b/Examples/Algorithms/Utilities/src/MakeMeasurementParticlesMap.cpp new file mode 100644 index 00000000000..b07573ddeda --- /dev/null +++ b/Examples/Algorithms/Utilities/src/MakeMeasurementParticlesMap.cpp @@ -0,0 +1,28 @@ +#include "ActsExamples/Utilities/MakeMeasurementParticlesMap.hpp" + +using namespace ActsExamples; + +ActsExamples::MakeMeasurementParticlesMap::MakeMeasurementParticlesMap( + Config cfg, Acts::Logging::Level lvl) + : IAlgorithm("MakeMeasurementParticlesMap", lvl), m_cfg(cfg) { + m_inputHitMap.initialize(m_cfg.inputMeasurementSimhitMap); + m_inputHits.initialize(m_cfg.inputSimHits); + m_outputParticleMap.initialize(m_cfg.outputMeasurementParticlesMap); +} + +ProcessCode ActsExamples::MakeMeasurementParticlesMap::execute( + const AlgorithmContext &ctx) const { + const auto hits = m_inputHits(ctx); + const auto hitMeasMap = m_inputHitMap(ctx); + + IndexMultimap outputMap; + + for (const auto &[measIdx, hitIdx] : hitMeasMap) { + const auto &hit = hits.nth(hitIdx); + outputMap.emplace(measIdx, hit->particleId()); + } + + m_outputParticleMap(ctx, std::move(outputMap)); + + return ProcessCode::SUCCESS; +} diff --git a/Examples/Framework/CMakeLists.txt b/Examples/Framework/CMakeLists.txt index 9ee5d487ddb..2752cd8998e 100644 --- a/Examples/Framework/CMakeLists.txt +++ b/Examples/Framework/CMakeLists.txt @@ -7,6 +7,7 @@ add_library( ActsExamplesFramework SHARED src/EventData/MeasurementCalibration.cpp src/EventData/ScalingCalibrator.cpp + src/EventData/GeometryContainers.cpp src/Framework/IAlgorithm.cpp src/Framework/SequenceElement.cpp src/Framework/WhiteBoard.cpp diff --git a/Examples/Framework/include/ActsExamples/EventData/Cluster.hpp b/Examples/Framework/include/ActsExamples/EventData/Cluster.hpp index 8fac374a4de..90565e304e1 100644 --- a/Examples/Framework/include/ActsExamples/EventData/Cluster.hpp +++ b/Examples/Framework/include/ActsExamples/EventData/Cluster.hpp @@ -8,7 +8,7 @@ #pragma once -#include "ActsFatras/Digitization/Channelizer.hpp" +#include "ActsFatras/Digitization/Segmentizer.hpp" #include @@ -16,7 +16,7 @@ namespace ActsExamples { /// Simple struct holding cluster information. struct Cluster { - using Cell = ActsFatras::Channelizer::ChannelSegment; + using Cell = ActsFatras::Segmentizer::ChannelSegment; size_t sizeLoc0 = 0; size_t sizeLoc1 = 0; std::vector channels; diff --git a/Examples/Framework/include/ActsExamples/EventData/GeometryContainers.hpp b/Examples/Framework/include/ActsExamples/EventData/GeometryContainers.hpp index e3007901bfc..fe01623c0ae 100644 --- a/Examples/Framework/include/ActsExamples/EventData/GeometryContainers.hpp +++ b/Examples/Framework/include/ActsExamples/EventData/GeometryContainers.hpp @@ -11,6 +11,7 @@ #include "Acts/EventData/SourceLink.hpp" #include "Acts/Geometry/GeometryIdentifier.hpp" #include "Acts/Surfaces/Surface.hpp" +#include "ActsExamples/EventData/Measurement.hpp" #include "ActsExamples/Utilities/GroupBy.hpp" #include "ActsExamples/Utilities/Range.hpp" @@ -56,6 +57,9 @@ struct GeometryIdGetter { -> decltype(thing.get().geometryId(), Acts::GeometryIdentifier()) { return thing.get().geometryId(); } + // support measurements (Implemented in cpp to avoid cyclic include) + Acts::GeometryIdentifier operator()( + const ActsExamples::Measurement& meas) const; }; struct CompareGeometryId { @@ -64,7 +68,7 @@ struct CompareGeometryId { // compare two elements using the automatic key extraction. template constexpr bool operator()(Left&& lhs, Right&& rhs) const { - return GeometryIdGetter()(lhs) < GeometryIdGetter()(rhs); + return GeometryIdGetter{}(lhs) < GeometryIdGetter{}(rhs); } }; diff --git a/Examples/Framework/src/EventData/GeometryContainers.cpp b/Examples/Framework/src/EventData/GeometryContainers.cpp new file mode 100644 index 00000000000..9bc195cbdd3 --- /dev/null +++ b/Examples/Framework/src/EventData/GeometryContainers.cpp @@ -0,0 +1,11 @@ +#include "ActsExamples/EventData/GeometryContainers.hpp" + +#include "ActsExamples/EventData/IndexSourceLink.hpp" + +Acts::GeometryIdentifier ActsExamples::detail::GeometryIdGetter::operator()( + const ActsExamples::Measurement& meas) const { + auto f = [](const auto& m) { + return m.sourceLink().template get().geometryId(); + }; + return std::visit(f, meas); +} diff --git a/Examples/Io/Csv/CMakeLists.txt b/Examples/Io/Csv/CMakeLists.txt index 0b5dd42d028..af70b4f0315 100644 --- a/Examples/Io/Csv/CMakeLists.txt +++ b/Examples/Io/Csv/CMakeLists.txt @@ -17,6 +17,7 @@ add_library( src/CsvTrackWriter.cpp src/CsvProtoTrackWriter.cpp src/CsvSpacePointWriter.cpp + src/CsvExaTrkXGraphWriter.cpp src/CsvBFieldWriter.cpp) target_include_directories( ActsExamplesIoCsv diff --git a/Examples/Io/Csv/include/ActsExamples/Io/Csv/CsvExaTrkXGraphWriter.hpp b/Examples/Io/Csv/include/ActsExamples/Io/Csv/CsvExaTrkXGraphWriter.hpp new file mode 100644 index 00000000000..9a420f90b23 --- /dev/null +++ b/Examples/Io/Csv/include/ActsExamples/Io/Csv/CsvExaTrkXGraphWriter.hpp @@ -0,0 +1,57 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2020 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/Utilities/Logger.hpp" +#include "ActsExamples/Framework/ProcessCode.hpp" +#include "ActsExamples/Framework/WriterT.hpp" +#include "ActsExamples/Utilities/Paths.hpp" + +#include +#include +#include + +namespace ActsExamples { +struct AlgorithmContext; + +class CsvExaTrkXGraphWriter final + : public WriterT, std::vector>> { + public: + struct Config { + /// Which simulated (truth) hits collection to use. + std::string inputGraph; + /// Where to place output files + std::string outputDir; + /// Output filename stem. + std::string outputStem = "exatrkx-graph"; + }; + + /// Construct the cluster writer. + /// + /// @param config is the configuration object + /// @param level is the logging level + CsvExaTrkXGraphWriter(const Config& config, Acts::Logging::Level level); + + /// Readonly access to the config + const Config& config() const { return m_cfg; } + + protected: + /// Type-specific write implementation. + /// + /// @param[in] ctx is the algorithm context + /// @param[in] simHits are the simhits to be written + ProcessCode writeT(const AlgorithmContext& ctx, + const std::pair, std::vector>& + graph) override; + + private: + Config m_cfg; +}; + +} // namespace ActsExamples diff --git a/Examples/Io/Csv/src/CsvExaTrkXGraphWriter.cpp b/Examples/Io/Csv/src/CsvExaTrkXGraphWriter.cpp new file mode 100644 index 00000000000..9ee9998b524 --- /dev/null +++ b/Examples/Io/Csv/src/CsvExaTrkXGraphWriter.cpp @@ -0,0 +1,56 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2020 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "ActsExamples/Io/Csv/CsvExaTrkXGraphWriter.hpp" + +#include "Acts/Definitions/Algebra.hpp" +#include "Acts/Definitions/Common.hpp" +#include "Acts/Definitions/Units.hpp" +#include "ActsExamples/Framework/AlgorithmContext.hpp" +#include "ActsExamples/Utilities/Paths.hpp" +#include "ActsFatras/EventData/Barcode.hpp" + +#include +#include + +#include +#include + +struct GraphData { + int64_t edge0; + int64_t edge1; + float weight; + DFE_NAMEDTUPLE(GraphData, edge0, edge1, weight); +}; + +ActsExamples::CsvExaTrkXGraphWriter::CsvExaTrkXGraphWriter( + const ActsExamples::CsvExaTrkXGraphWriter::Config& config, + Acts::Logging::Level level) + : WriterT(config.inputGraph, "CsvExaTrkXGraphWriter", level), + m_cfg(config) {} + +ActsExamples::ProcessCode ActsExamples::CsvExaTrkXGraphWriter::writeT( + const ActsExamples::AlgorithmContext& ctx, + const std::pair, std::vector>& graph) { + std::string path = perEventFilepath( + m_cfg.outputDir, m_cfg.outputStem + ".csv", ctx.eventNumber); + + dfe::NamedTupleCsvWriter writer(path); + + const auto& [edges, weights] = graph; + + for (auto i = 0ul; i < weights.size(); ++i) { + GraphData edge; + edge.edge0 = edges[2 * i]; + edge.edge1 = edges[2 * i + 1]; + edge.weight = weights[i]; + writer.append(edge); + } + + return ActsExamples::ProcessCode::SUCCESS; +} diff --git a/Examples/Io/Csv/src/CsvMeasurementReader.cpp b/Examples/Io/Csv/src/CsvMeasurementReader.cpp index 67b905e81d7..bc09e90cd9a 100644 --- a/Examples/Io/Csv/src/CsvMeasurementReader.cpp +++ b/Examples/Io/Csv/src/CsvMeasurementReader.cpp @@ -145,10 +145,10 @@ ActsExamples::ClusterContainer makeClusters( for (auto it = begin; it != end; ++it) { const auto& cellData = it->second; - ActsFatras::Channelizer::Segment2D dummySegment = {Acts::Vector2::Zero(), + ActsFatras::Segmentizer::Segment2D dummySegment = {Acts::Vector2::Zero(), Acts::Vector2::Zero()}; - ActsFatras::Channelizer::Bin2D bin{ + ActsFatras::Segmentizer::Bin2D bin{ static_cast(cellData.channel0), static_cast(cellData.channel1)}; diff --git a/Examples/Io/Csv/src/CsvParticleReader.cpp b/Examples/Io/Csv/src/CsvParticleReader.cpp index 6c4dd306196..4a52ddd919d 100644 --- a/Examples/Io/Csv/src/CsvParticleReader.cpp +++ b/Examples/Io/Csv/src/CsvParticleReader.cpp @@ -59,6 +59,8 @@ ActsExamples::ProcessCode ActsExamples::CsvParticleReader::read( auto path = perEventFilepath(m_cfg.inputDir, m_cfg.inputStem + ".csv", ctx.eventNumber); + ACTS_DEBUG("Read '" << path << "'"); + // vt and m are an optional columns dfe::NamedTupleCsvReader reader(path, {"vt", "m"}); ParticleData data; diff --git a/Examples/Io/Csv/src/CsvSimHitReader.cpp b/Examples/Io/Csv/src/CsvSimHitReader.cpp index 9db7cac3aa6..d75cb570d42 100644 --- a/Examples/Io/Csv/src/CsvSimHitReader.cpp +++ b/Examples/Io/Csv/src/CsvSimHitReader.cpp @@ -54,6 +54,7 @@ ActsExamples::ProcessCode ActsExamples::CsvSimHitReader::read( const ActsExamples::AlgorithmContext& ctx) { auto path = perEventFilepath(m_cfg.inputDir, m_cfg.inputStem + ".csv", ctx.eventNumber); + ACTS_DEBUG("Read '" << path << "'"); dfe::NamedTupleCsvReader reader(path); diff --git a/Examples/Io/Performance/ActsExamples/Io/Performance/CKFPerformanceWriter.cpp b/Examples/Io/Performance/ActsExamples/Io/Performance/CKFPerformanceWriter.cpp index 548ee6dbd69..4a0fc4028fb 100644 --- a/Examples/Io/Performance/ActsExamples/Io/Performance/CKFPerformanceWriter.cpp +++ b/Examples/Io/Performance/ActsExamples/Io/Performance/CKFPerformanceWriter.cpp @@ -19,6 +19,7 @@ #include #include +#include #include #include #include @@ -68,6 +69,8 @@ ActsExamples::CKFPerformanceWriter::CKFPerformanceWriter( m_fakeRatePlotTool.book(m_fakeRatePlotCache); m_duplicationPlotTool.book(m_duplicationPlotCache); m_trackSummaryPlotTool.book(m_trackSummaryPlotCache); + + m_particlesMatchedStream << "event,particle_id,matched\n"; } ActsExamples::CKFPerformanceWriter::~CKFPerformanceWriter() { @@ -130,6 +133,12 @@ ActsExamples::ProcessCode ActsExamples::CKFPerformanceWriter::finalize() { write_float(duplicationRate_particle, "duplicaterate_particles"); ACTS_INFO("Wrote performance plots to '" << m_outputFile->GetPath() << "'"); } + + const auto csvFileName = + m_cfg.filePath.substr(0, m_cfg.filePath.size() - 5) + ".csv"; + std::ofstream outputFileParticlesMatched(csvFileName); + outputFileParticlesMatched << m_particlesMatchedStream.str(); + return ProcessCode::SUCCESS; } @@ -304,5 +313,11 @@ ActsExamples::ProcessCode ActsExamples::CKFPerformanceWriter::writeT( m_nTotalParticles += 1; } // end all truth particles + for (const auto& p : particles) { + bool m = (matched.find(p.particleId()) != matched.end()); + m_particlesMatchedStream << ctx.eventNumber << "," << p.particleId().value() + << "," << m << "\n"; + } + return ProcessCode::SUCCESS; } diff --git a/Examples/Io/Performance/ActsExamples/Io/Performance/CKFPerformanceWriter.hpp b/Examples/Io/Performance/ActsExamples/Io/Performance/CKFPerformanceWriter.hpp index 1ae1ae39c99..a6f1c2470c9 100644 --- a/Examples/Io/Performance/ActsExamples/Io/Performance/CKFPerformanceWriter.hpp +++ b/Examples/Io/Performance/ActsExamples/Io/Performance/CKFPerformanceWriter.hpp @@ -121,6 +121,8 @@ class CKFPerformanceWriter final : public WriterT { ReadDataHandle m_inputParticles{this, "InputParticles"}; ReadDataHandle m_inputMeasurementParticlesMap{ this, "InputMeasurementParticlesMap"}; + + std::stringstream m_particlesMatchedStream; }; } // namespace ActsExamples diff --git a/Examples/Io/Performance/ActsExamples/Io/Performance/SeedingPerformanceWriter.cpp b/Examples/Io/Performance/ActsExamples/Io/Performance/SeedingPerformanceWriter.cpp index d63037308ff..fce18c458df 100644 --- a/Examples/Io/Performance/ActsExamples/Io/Performance/SeedingPerformanceWriter.cpp +++ b/Examples/Io/Performance/ActsExamples/Io/Performance/SeedingPerformanceWriter.cpp @@ -22,6 +22,8 @@ #include #include +#include +#include namespace ActsExamples { struct AlgorithmContext; @@ -95,10 +97,22 @@ ActsExamples::ProcessCode ActsExamples::SeedingPerformanceWriter::finalize() { "/ nMatchedParticles) = " << aveNDuplicatedSeeds); + auto write_float = [&](float f, const char* name) { + TVectorF v(1); + v[0] = f; + m_outputFile->WriteObject(&v, name); + }; + if (m_outputFile != nullptr) { m_outputFile->cd(); m_effPlotTool.write(m_effPlotCache); m_duplicationPlotTool.write(m_duplicationPlotCache); + write_float(m_nTotalSeeds, "total_seeds"); + write_float(totalSeedPurity, "total_seed_purity"); + write_float(eff, "seed_efficiency"); + write_float(fakeRate, "seed_fakerate"); + write_float(duplicationRate, "seed_duplicationrate"); + write_float(aveNDuplicatedSeeds, "avg_duplicate_seeds"); ACTS_INFO("Wrote performance plots to '" << m_outputFile->GetPath() << "'"); } return ProcessCode::SUCCESS; diff --git a/Examples/Python/python/acts/examples/__init__.py b/Examples/Python/python/acts/examples/__init__.py index 0c37423ce90..91e8061fe81 100644 --- a/Examples/Python/python/acts/examples/__init__.py +++ b/Examples/Python/python/acts/examples/__init__.py @@ -37,6 +37,7 @@ def ConcretePropagator(propagator): "OnnxMetricLearning", "TorchEdgeClassifier", "OnnxEdgeClassifier", + "BoostTrackBuilding", ]: if hasattr(ActsPythonBindings._examples, module): _patchKwargsConstructor(getattr(ActsPythonBindings._examples, module)) @@ -527,9 +528,15 @@ def _getAutoFpeMasks(cls) -> List[FpeMask]: @classmethod def _printFpeSummary(cls, masks: List[FpeMask]): + disableEnvironVar = "ACTS_SEQUENCER_DISABLE_FPE_MASK_PRINTING" + if disableEnvironVar in os.environ: + return + if len(masks) == 0: return + print(f"Note: Disable FPE mask printing by defining {disableEnvironVar}") + # Try to make a nice summary with rich, or fallback to a plain text one try: import rich diff --git a/Examples/Python/python/acts/examples/reconstruction.py b/Examples/Python/python/acts/examples/reconstruction.py index 36b311a6781..84ac623e1f3 100644 --- a/Examples/Python/python/acts/examples/reconstruction.py +++ b/Examples/Python/python/acts/examples/reconstruction.py @@ -1359,11 +1359,13 @@ def addExaTrkX( filterConfig = { "level": customLogLevel(), + "numFeatures": 3, "cut": 0.01, } gnnConfig = { "level": customLogLevel(), + "numFeatures": 3, "cut": 0.5, } diff --git a/Examples/Python/src/ExaTrkXTrackFinding.cpp b/Examples/Python/src/ExaTrkXTrackFinding.cpp index ad27045da03..d112943e044 100644 --- a/Examples/Python/src/ExaTrkXTrackFinding.cpp +++ b/Examples/Python/src/ExaTrkXTrackFinding.cpp @@ -18,6 +18,7 @@ #include "ActsExamples/TrackFindingExaTrkX/PrototracksToParameters.hpp" #include "ActsExamples/TrackFindingExaTrkX/TrackFindingAlgorithmExaTrkX.hpp" #include "ActsExamples/TrackFindingExaTrkX/TrackFindingFromPrototrackAlgorithm.hpp" +#include "ActsExamples/TrackFindingExaTrkX/ProtoTrackEffPurPrinter.hpp" #include @@ -70,6 +71,7 @@ void addExaTrkXTrackFinding(Context &ctx) { ACTS_PYTHON_MEMBER(embeddingDim); ACTS_PYTHON_MEMBER(rVal); ACTS_PYTHON_MEMBER(knnVal); + ACTS_PYTHON_MEMBER(shuffleDirections); ACTS_PYTHON_STRUCT_END(); } { @@ -103,8 +105,7 @@ void addExaTrkXTrackFinding(Context &ctx) { .def(py::init([](Logging::Level lvl) { return std::make_shared( getDefaultLogger("EdgeClassifier", lvl)); - }), - py::arg("level")); + }), py::arg("level")); } #endif @@ -169,9 +170,9 @@ void addExaTrkXTrackFinding(Context &ctx) { ActsExamples::TrackFindingAlgorithmExaTrkX, mex, "TrackFindingAlgorithmExaTrkX", inputSpacePoints, inputSimHits, inputParticles, inputClusters, inputMeasurementSimhitsMap, - outputProtoTracks, graphConstructor, edgeClassifiers, trackBuilder, - rScale, phiScale, zScale, cellCountScale, cellSumScale, clusterXScale, - clusterYScale, targetMinHits, targetMinPT); + outputProtoTracks, outputGraph, graphConstructor, edgeClassifiers, + trackBuilder, rScale, phiScale, zScale, cellCountScale, cellSumScale, + clusterXScale, clusterYScale, targetMinHits, targetMinPT); { auto cls = @@ -212,6 +213,10 @@ void addExaTrkXTrackFinding(Context &ctx) { py::arg("timing") = nullptr); } + ACTS_PYTHON_DECLARE_ALGORITHM(ActsExamples::ProtoTrackEffPurPrinter, mex, + "ProtoTrackEffPurPrinter", testProtoTracks, + refProtoTracks); + ACTS_PYTHON_DECLARE_ALGORITHM( ActsExamples::PrototracksToParameters, mex, "PrototracksToParameters", inputProtoTracks, inputSpacePoints, outputSeeds, outputParameters, diff --git a/Examples/Python/src/ModuleEntry.cpp b/Examples/Python/src/ModuleEntry.cpp index 0c8f3f7ee41..f269a153b8a 100644 --- a/Examples/Python/src/ModuleEntry.cpp +++ b/Examples/Python/src/ModuleEntry.cpp @@ -74,7 +74,6 @@ void addHepMC3(Context& ctx); void addExaTrkXTrackFinding(Context& ctx); void addEDM4hep(Context& ctx); void addSvg(Context& ctx); -void addObj(Context& ctx); void addOnnx(Context& ctx); void addOnnxMlpack(Context& ctx); void addOnnxNeuralCalibrator(Context& ctx); @@ -124,7 +123,6 @@ PYBIND11_MODULE(ActsPythonBindings, m) { addHepMC3(ctx); addExaTrkXTrackFinding(ctx); addEDM4hep(ctx); - addObj(ctx); addSvg(ctx); addOnnx(ctx); addOnnxMlpack(ctx); diff --git a/Examples/Python/src/Output.cpp b/Examples/Python/src/Output.cpp index e5dde8abcea..a544773c0fc 100644 --- a/Examples/Python/src/Output.cpp +++ b/Examples/Python/src/Output.cpp @@ -14,6 +14,7 @@ #include "ActsExamples/Digitization/DigitizationConfig.hpp" #include "ActsExamples/Framework/ProcessCode.hpp" #include "ActsExamples/Io/Csv/CsvBFieldWriter.hpp" +#include "ActsExamples/Io/Csv/CsvExaTrkXGraphWriter.hpp" #include "ActsExamples/Io/Csv/CsvMeasurementWriter.hpp" #include "ActsExamples/Io/Csv/CsvParticleWriter.hpp" #include "ActsExamples/Io/Csv/CsvPlanarClusterWriter.hpp" @@ -373,7 +374,7 @@ void addOutput(Context& ctx) { inputParticles, inputMeasurementParticlesMap, filePath, fileMode, effPlotToolConfig, fakeRatePlotToolConfig, duplicationPlotToolConfig, - trackSummaryPlotToolConfig, duplicatedPredictor); + trackSummaryPlotToolConfig, duplicatedPredictor, truthMatchProbMin, doubleMatching); ACTS_PYTHON_DECLARE_WRITER( ActsExamples::RootNuclearInteractionParametersWriter, mex, @@ -405,5 +406,9 @@ void addOutput(Context& ctx) { register_csv_bfield_writer_binding(w); register_csv_bfield_writer_binding(w); } + + ACTS_PYTHON_DECLARE_WRITER(ActsExamples::CsvExaTrkXGraphWriter, mex, + "CsvExaTrkXGraphWriter", inputGraph, outputDir, + outputStem); } } // namespace Acts::Python diff --git a/Examples/Python/src/TrackFinding.cpp b/Examples/Python/src/TrackFinding.cpp index f26cd169db9..8070f611a37 100644 --- a/Examples/Python/src/TrackFinding.cpp +++ b/Examples/Python/src/TrackFinding.cpp @@ -32,6 +32,14 @@ #include "ActsExamples/Utilities/TracksToParameters.hpp" #include "ActsExamples/Utilities/TracksToTrajectories.hpp" #include "ActsExamples/Utilities/TrajectoriesToPrototracks.hpp" +#include "ActsExamples/Utilities/MeasurementMapSelector.hpp" + +//////////////////////////// +// For GNN+CKF Experiment // +//////////////////////////// +// #include "ActsExamples/TrackFindingX/ParameterFromTrajectoryAlgorithm.hpp" +// #include "ActsExamples/TrackFindingX/SourceLinkSelectorAlgorithm.hpp" +#include "ActsExamples/Utilities/MakeMeasurementParticlesMap.hpp" #include #include @@ -388,6 +396,27 @@ void addTrackFinding(Context& ctx) { ActsExamples::MeasurementMapSelector, mex, "MeasurementMapSelector", inputMeasurementParticleMap, inputSourceLinks, outputMeasurementParticleMap, geometrySelection); + + //////////////////////////// + // For GNN+CKF Experiment // + //////////////////////////// + { + // ACTS_PYTHON_DECLARE_ALGORITHM(ActsExamples::SourceLinkSelectorAlgorithm, + // mex, + // "SourceLinkSelectorAlgorithm", + // inputSourceLinks, outputSourceLinks, + // geometrySelection); + // + // ACTS_PYTHON_DECLARE_ALGORITHM( + // ActsExamples::ParameterFromTrajectoryAlgorithm, mex, + // "ParameterFromTrajectoryAlgorithm", inputTrajectories, + // outputParamters); + + ACTS_PYTHON_DECLARE_ALGORITHM(ActsExamples::MakeMeasurementParticlesMap, + mex, "MakeMeasurementParticlesMap", + inputMeasurementSimhitMap, inputSimHits, + outputMeasurementParticlesMap); + } } } // namespace Acts::Python diff --git a/Examples/Python/src/TruthTracking.cpp b/Examples/Python/src/TruthTracking.cpp index a49a2e972d7..dec7875592d 100644 --- a/Examples/Python/src/TruthTracking.cpp +++ b/Examples/Python/src/TruthTracking.cpp @@ -17,6 +17,7 @@ #include "ActsExamples/TruthTracking/TruthSeedingAlgorithm.hpp" #include "ActsExamples/TruthTracking/TruthTrackFinder.hpp" #include "ActsExamples/TruthTracking/TruthVertexFinder.hpp" +#include "ActsExamples/Utilities/HitSelector.hpp" #include "ActsExamples/Utilities/Range.hpp" #include @@ -42,7 +43,7 @@ void addTruthTracking(Context& ctx) { ACTS_PYTHON_DECLARE_ALGORITHM( ActsExamples::TruthTrackFinder, mex, "TruthTrackFinder", inputParticles, - inputMeasurementParticlesMap, outputProtoTracks); + inputMeasurementParticlesMap, outputProtoTracks, minHits); { using Alg = ActsExamples::TruthSeedSelector; @@ -195,6 +196,9 @@ void addTruthTracking(Context& ctx) { ActsExamples::TruthSeedingAlgorithm, mex, "TruthSeedingAlgorithm", inputParticles, inputMeasurementParticlesMap, inputSpacePoints, outputParticles, outputSeeds, outputProtoTracks, deltaRMin, deltaRMax); + + ACTS_PYTHON_DECLARE_ALGORITHM(ActsExamples::HitSelector, mex, "HitSelector", + inputHits, outputHits, maxTime); } } // namespace Acts::Python diff --git a/Examples/Scripts/Python/exatrkx.py b/Examples/Scripts/Python/exatrkx.py index 16cc606de20..590db2def79 100755 --- a/Examples/Scripts/Python/exatrkx.py +++ b/Examples/Scripts/Python/exatrkx.py @@ -79,6 +79,7 @@ modelDir, outputDir, backend=backend, + logLevel=acts.logging.VERBOSE ) s.run() diff --git a/Examples/Scripts/Python/exatrkx_plain_pipeline.py b/Examples/Scripts/Python/exatrkx_plain_pipeline.py new file mode 100755 index 00000000000..8f712440661 --- /dev/null +++ b/Examples/Scripts/Python/exatrkx_plain_pipeline.py @@ -0,0 +1,115 @@ +#!/usr/bin/env python3 + +from multiprocessing import Process, Queue +from pathlib import Path +import argparse +import os + +import numpy as np + + +def load_pytorch_graph(queue, filename): + import torch + + graph = torch.load(filename) + + graphdict = {} + for k in graph.keys: + try: + graphdict[k] = graph[k].detach().cpu().numpy() + except: + continue + + queue.put(graphdict) + + +def run_pipeline(data, torchscript_dir): + import acts + import acts.examples + + # Make input data + scaling = ScalingDict( + { + "r": 1000.0, + "phi": 3.14159, + "z": 3000.0, + } + ) + + input_tensor = np.vstack( + [ + data[k] / scaling[k] + for k in ["r", "phi", "z", "cell_count", "cell_val", "lx", "ly"] + ] + ).T + + # Make metric hook + hook = acts.examples.TorchTruthGraphMetricsHook( + data["track_edges"].T.flatten(), acts.logging.DEBUG + ) + + # Make stages + emb = acts.examples.TorchMetricLearning( + acts.logging.VERBOSE, + embeddingDim=8, + knnVal=100, + numFeatures=7, + rVal=0.2, + modelPath=torchscript_dir / "embedding.pt", + ) + + flt = acts.examples.TorchEdgeClassifier( + acts.logging.VERBOSE, + numFeatures=3, + undirected=False, + modelPath=torchscript_dir / "filter.pt", + ) + + gnn = acts.examples.TorchEdgeClassifier( + acts.logging.VERBOSE, + undirected=True, + numFeatures=3, + modelPath=torchscript_dir / "gnn.pt", + ) + + trk = acts.examples.BoostTrackBuilding(acts.logging.VERBOSE) + + # Run pipeline + pipeline = acts.examples.Pipeline(emb, [flt,gnn], trk, acts.logging.VERBOSE) + + return pipeline.run( + input_tensor.flatten().tolist(), np.arange(input_tensor.shape[0]).tolist(), hook + ) + + +class ScalingDict(dict): + def __init__(self, *args, **kwargs): + super().__init__(*args, **kwargs) + + def __missing__(self, key): + return 1.0 + + +def main(): + + parser = argparse.ArgumentParser() + parser.add_argument("torchscript_dir", type=str) + parser.add_argument("pyg_file", type=str) + args = vars(parser.parse_args()) + + assert os.path.exists(args["torchscript_dir"]) + assert os.path.exists(args["pyg_file"]) + + # Load *.pyg in different process to avoid conflicting torch versions for python and C++ + q = Queue() + p = Process(target=load_pytorch_graph, args=(q, args["pyg_file"])) + p.start() + graph = q.get() + p.join() + + run_pipeline(graph, Path(args["torchscript_dir"])) + + + +if __name__ == "__main__": + main() diff --git a/Fatras/CMakeLists.txt b/Fatras/CMakeLists.txt index f05fb13ccc7..0edc77daf99 100644 --- a/Fatras/CMakeLists.txt +++ b/Fatras/CMakeLists.txt @@ -1,6 +1,6 @@ add_library( ActsFatras SHARED - src/Digitization/Channelizer.cpp + src/Digitization/Segmentizer.cpp src/Digitization/DigitizationError.cpp src/Digitization/PlanarSurfaceMask.cpp src/Digitization/PlanarSurfaceDrift.cpp diff --git a/Fatras/include/ActsFatras/Digitization/Channelizer.hpp b/Fatras/include/ActsFatras/Digitization/Channelizer.hpp index c40cfc1d5c7..199401bc550 100644 --- a/Fatras/include/ActsFatras/Digitization/Channelizer.hpp +++ b/Fatras/include/ActsFatras/Digitization/Channelizer.hpp @@ -8,97 +8,61 @@ #pragma once -#include "Acts/Definitions/Algebra.hpp" -#include "Acts/Geometry/GeometryContext.hpp" +#include "ActsFatras/Digitization/PlanarSurfaceDrift.hpp" +#include "ActsFatras/Digitization/PlanarSurfaceMask.hpp" +#include "ActsFatras/Digitization/Segmentizer.hpp" +#include "ActsFatras/EventData/Hit.hpp" -#include -#include -#include - -namespace Acts { -class BinUtility; -class Surface; -} // namespace Acts +#include namespace ActsFatras { -/// The Channelizer splits a surface segment, i.e. after projection -/// onto the readout surface into channel segments. -/// -struct Channelizer { - /// Shorthand for a 2D segment - using Segment2D = std::array; - /// Shorthand for a 2D bin - using Bin2D = std::array; - /// shorthand for a 2D bin delta - using BinDelta2D = std::array; +/// @brief Class that ties the digitization modules together and produces the channels +class Channelizer { + PlanarSurfaceDrift m_surfaceDrift; + PlanarSurfaceMask m_surfaceMask; + Segmentizer m_segmentizer; - /// Nested struct for stepping from one channel to the next. - struct ChannelStep { - /// This is the delta to the last step in bins - BinDelta2D delta = {0, 0}; - /// The intersection with the channel boundary - Acts::Vector2 intersect; - /// The patlength from the start - double path = 0.; + public: + /// Do the geometric channelizing + /// + /// @param geoCfg is the geometric digitization configuration + /// @param hit the Simultated hit + /// @param surface the Surface on which this is supposed to happen + /// @param gctx the Geometry context + /// + /// @return the list of channels + Acts::Result> channelize( + const Hit& hit, const Acts::Surface& surface, + const Acts::GeometryContext& gctx, const Acts::Vector3& driftDir, + const Acts::BinUtility& segmentation, double thickness) const { + auto driftedSegment = m_surfaceDrift.toReadout( + gctx, surface, thickness, hit.position(), hit.direction(), driftDir); - /// Constructor with arguments for a ChannelStep. - /// - /// @param delta_ The bin delta for this step - /// @param intersect_ The intersect with the channel boundary - /// @param start The start of the surface segment, for path from origin - ChannelStep(BinDelta2D delta_, Acts::Vector2 intersect_, - const Acts::Vector2& start) - : delta(delta_), - intersect(std::move(intersect_)), - path((intersect - start).norm()) {} + auto maskedSegmentRes = m_surfaceMask.apply(surface, driftedSegment); - /// Smaller operator for sorting the ChannelStep objects. - /// - /// @param cstep The other ChannelStep to be compared - /// - /// The ChannelStep objects can be compared with its path distance - /// from the start (surface segment origin) - bool operator<(const ChannelStep& cstep) const { return path < cstep.path; } - }; + if (!maskedSegmentRes.ok()) { + return maskedSegmentRes.error(); + } - /// Nested struct for representing channel steps. - struct ChannelSegment { - /// The bin of this segment - Bin2D bin = {0, 0}; - /// The segment start, end points - Segment2D path2D; - /// The (clipped) value (uncorrected: path length) - double activation = 0.; + // Now Channelize + auto segments = + m_segmentizer.segments(gctx, surface, segmentation, *maskedSegmentRes); - /// Constructor with arguments - /// - /// @param bin_ The bin corresponding to this step - /// @param path2D_ The start/end 2D position of the segement - /// @param activation_ The segment activation (clean: length) for this bin - ChannelSegment(Bin2D bin_, Segment2D path2D_, double activation_) - : bin(bin_), path2D(std::move(path2D_)), activation(activation_) {} - }; + // Go from 2D-path to 3D-path by applying thickness + const auto path2D = std::accumulate( + segments.begin(), segments.end(), 0.0, + [](double sum, const auto& seg) { return sum + seg.activation; }); - /// Divide the surface segment into channel segments. - /// - /// @note Channelizing is done in cartesian coordinates (start/end) - /// @note The start and end cartesian vector is supposed to be inside - /// the surface bounds (pre-run through the SurfaceMasker) - /// @note The segmentation has to be 2-dimensional, even if the - /// actual readout is 1-dimensional, in latter case one bin in the - /// second coordinate direction is required. - /// - /// @param geoCtx The geometry context for the localToGlobal, etc. - /// @param surface The surface for the channelizing - /// @param segmentation The segmentation for the channelizing - /// @param segment The surface segment (cartesian coordinates) - /// - /// @return a vector of ChannelSegment objects - std::vector segments(const Acts::GeometryContext& geoCtx, - const Acts::Surface& surface, - const Acts::BinUtility& segmentation, - const Segment2D& segment) const; + for (auto& seg : segments) { + auto r = path2D != 0.0 ? (seg.activation / path2D) : 1.0; + auto segThickness = r * thickness; + + seg.activation = std::hypot(segThickness, seg.activation); + } + + return segments; + } }; } // namespace ActsFatras diff --git a/Fatras/include/ActsFatras/Digitization/PlanarSurfaceMask.hpp b/Fatras/include/ActsFatras/Digitization/PlanarSurfaceMask.hpp index 6474ddc086f..1202c51a568 100644 --- a/Fatras/include/ActsFatras/Digitization/PlanarSurfaceMask.hpp +++ b/Fatras/include/ActsFatras/Digitization/PlanarSurfaceMask.hpp @@ -6,6 +6,8 @@ // License, v. 2.0. If a copy of the MPL was not distributed with this // file, You can obtain one at http://mozilla.org/MPL/2.0/. +#pragma once + #include "Acts/Definitions/Algebra.hpp" #include "Acts/Geometry/GeometryContext.hpp" #include "Acts/Surfaces/detail/IntersectionHelper2D.hpp" diff --git a/Fatras/include/ActsFatras/Digitization/Segmentizer.hpp b/Fatras/include/ActsFatras/Digitization/Segmentizer.hpp new file mode 100644 index 00000000000..e6e8ac6327d --- /dev/null +++ b/Fatras/include/ActsFatras/Digitization/Segmentizer.hpp @@ -0,0 +1,104 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2020 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/Definitions/Algebra.hpp" +#include "Acts/Geometry/GeometryContext.hpp" + +#include +#include +#include + +namespace Acts { +class BinUtility; +class Surface; +} // namespace Acts + +namespace ActsFatras { + +/// The Segmentizer splits a surface segment, i.e. after projection +/// onto the readout surface into channel segments. +/// +struct Segmentizer { + /// Shorthand for a 2D segment + using Segment2D = std::array; + /// Shorthand for a 2D bin + using Bin2D = std::array; + /// shorthand for a 2D bin delta + using BinDelta2D = std::array; + + /// Nested struct for stepping from one channel to the next. + struct ChannelStep { + /// This is the delta to the last step in bins + BinDelta2D delta = {0, 0}; + /// The intersection with the channel boundary + Acts::Vector2 intersect; + /// The patlength from the start + double path = 0.; + + /// Constructor with arguments for a ChannelStep. + /// + /// @param delta_ The bin delta for this step + /// @param intersect_ The intersect with the channel boundary + /// @param start The start of the surface segment, for path from origin + ChannelStep(BinDelta2D delta_, Acts::Vector2 intersect_, + const Acts::Vector2& start) + : delta(delta_), + intersect(std::move(intersect_)), + path((intersect - start).norm()) {} + + /// Smaller operator for sorting the ChannelStep objects. + /// + /// @param cstep The other ChannelStep to be compared + /// + /// The ChannelStep objects can be compared with its path distance + /// from the start (surface segment origin) + bool operator<(const ChannelStep& cstep) const { return path < cstep.path; } + }; + + /// Nested struct for representing channel steps. + struct ChannelSegment { + /// The bin of this segment + Bin2D bin = {0, 0}; + /// The segment start, end points + Segment2D path2D; + /// The (clipped) value (uncorrected: path length) + double activation = 0.; + + /// Constructor with arguments + /// + /// @param bin_ The bin corresponding to this step + /// @param path2D_ The start/end 2D position of the segement + /// @param activation_ The segment activation (clean: length) for this bin + ChannelSegment(Bin2D bin_, Segment2D path2D_, double activation_) + : bin(bin_), path2D(std::move(path2D_)), activation(activation_) {} + }; + + /// Divide the surface segment into channel segments. + /// + /// @note Channelizing is done in cartesian coordinates (start/end) + /// @note The start and end cartesian vector is supposed to be inside + /// the surface bounds (pre-run through the SurfaceMasker) + /// @note The segmentation has to be 2-dimensional, even if the + /// actual readout is 1-dimensional, in latter case one bin in the + /// second coordinate direction is required. + /// + /// @param geoCtx The geometry context for the localToGlobal, etc. + /// @param surface The surface for the channelizing + /// @param segmentation The segmentation for the channelizing + /// @param segment The surface segment (cartesian coordinates) + /// + /// @return a vector of ChannelSegment objects + std::vector segments(const Acts::GeometryContext& geoCtx, + const Acts::Surface& surface, + const Acts::BinUtility& segmentation, + const Segment2D& segment) const; +}; + +} // namespace ActsFatras diff --git a/Fatras/src/Digitization/Channelizer.cpp b/Fatras/src/Digitization/Channelizer.cpp deleted file mode 100644 index fd541f3504e..00000000000 --- a/Fatras/src/Digitization/Channelizer.cpp +++ /dev/null @@ -1,164 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2020 CERN for the benefit of the Acts project -// -// This Source Code Form is subject to the terms of the Mozilla Public -// License, v. 2.0. If a copy of the MPL was not distributed with this -// file, You can obtain one at http://mozilla.org/MPL/2.0/. - -#include "ActsFatras/Digitization/Channelizer.hpp" - -#include "Acts/Surfaces/Surface.hpp" -#include "Acts/Surfaces/detail/IntersectionHelper2D.hpp" -#include "Acts/Utilities/BinUtility.hpp" -#include "Acts/Utilities/BinningType.hpp" -#include "Acts/Utilities/Helpers.hpp" -#include "Acts/Utilities/Intersection.hpp" - -#include -#include -#include - -std::vector -ActsFatras::Channelizer::segments(const Acts::GeometryContext& geoCtx, - const Acts::Surface& surface, - const Acts::BinUtility& segmentation, - const Segment2D& segment) const { - // Return if the segmentation is not two-dimensional - // (strips need to have one bin along the strip) - if (segmentation.dimensions() != 2) { - return {}; - } - - // Start and end point - const auto& start = segment[0]; - const auto& end = segment[1]; - - // Full path length - the full channel - auto segment2d = (end - start); - std::vector cSteps; - Bin2D bstart = {0, 0}; - Bin2D bend = {0, 0}; - - if (surface.type() == Acts::Surface::SurfaceType::Plane) { - // Get the segmentation and convert it to lines & arcs - bstart = {static_cast(segmentation.bin(start, 0)), - static_cast(segmentation.bin(start, 1))}; - bend = {static_cast(segmentation.bin(end, 0)), - static_cast(segmentation.bin(end, 1))}; - // Fast single channel exit - if (bstart == bend) { - return {ChannelSegment(bstart, {start, end}, segment2d.norm())}; - } - // The lines channel segment lines along x - if (bstart[0] != bend[0]) { - double k = segment2d.y() / segment2d.x(); - double d = start.y() - k * start.x(); - - const auto& xboundaries = segmentation.binningData()[0].boundaries(); - std::vector xbbounds = { - xboundaries.begin() + std::min(bstart[0], bend[0]) + 1, - xboundaries.begin() + std::max(bstart[0], bend[0]) + 1}; - for (const auto x : xbbounds) { - cSteps.push_back(ChannelStep{ - {(bstart[0] < bend[0] ? 1 : -1), 0}, {x, k * x + d}, start}); - } - } - // The lines channel segment lines along y - if (bstart[1] != bend[1]) { - double k = segment2d.x() / segment2d.y(); - double d = start.x() - k * start.y(); - const auto& yboundaries = segmentation.binningData()[1].boundaries(); - std::vector ybbounds = { - yboundaries.begin() + std::min(bstart[1], bend[1]) + 1, - yboundaries.begin() + std::max(bstart[1], bend[1]) + 1}; - for (const auto y : ybbounds) { - cSteps.push_back(ChannelStep{ - {0, (bstart[1] < bend[1] ? 1 : -1)}, {k * y + d, y}, start}); - } - } - - } else if (surface.type() == Acts::Surface::SurfaceType::Disc) { - Acts::Vector2 pstart(Acts::VectorHelpers::perp(start), - Acts::VectorHelpers::phi(start)); - Acts::Vector2 pend(Acts::VectorHelpers::perp(end), - Acts::VectorHelpers::phi(end)); - - // Get the segmentation and convert it to lines & arcs - bstart = {static_cast(segmentation.bin(pstart, 0)), - static_cast(segmentation.bin(pstart, 1))}; - bend = {static_cast(segmentation.bin(pend, 0)), - static_cast(segmentation.bin(pend, 1))}; - - // Fast single channel exit - if (bstart == bend) { - return {ChannelSegment(bstart, {start, end}, segment2d.norm())}; - } - - double phistart = pstart[1]; - double phiend = pend[1]; - - // The radial boundaries - if (bstart[0] != bend[0]) { - const auto& rboundaries = segmentation.binningData()[0].boundaries(); - std::vector rbbounds = { - rboundaries.begin() + std::min(bstart[0], bend[0]) + 1, - rboundaries.begin() + std::max(bstart[0], bend[0]) + 1}; - for (const auto& r : rbbounds) { - auto radIntersection = - Acts::detail::IntersectionHelper2D::intersectCircleSegment( - r, std::min(phistart, phiend), std::max(phistart, phiend), - start, (end - start).normalized()); - cSteps.push_back(ChannelStep{{(bstart[0] < bend[0] ? 1 : -1), 0}, - radIntersection.position(), - start}); - } - } - // The phi boundaries - if (bstart[1] != bend[1]) { - double referenceR = surface.binningPositionValue(geoCtx, Acts::binR); - Acts::Vector2 origin = {0., 0.}; - const auto& phiboundaries = segmentation.binningData()[1].boundaries(); - std::vector phibbounds = { - phiboundaries.begin() + std::min(bstart[1], bend[1]) + 1, - phiboundaries.begin() + std::max(bstart[1], bend[1]) + 1}; - - for (const auto& phi : phibbounds) { - Acts::Vector2 philine(referenceR * std::cos(phi), - referenceR * std::sin(phi)); - auto phiIntersection = - Acts::detail::IntersectionHelper2D::intersectSegment( - origin, philine, start, (end - start).normalized()); - cSteps.push_back(ChannelStep{{0, (bstart[1] < bend[1] ? 1 : -1)}, - phiIntersection.position(), - start}); - } - } - } - - // Register the last step if successful - if (!cSteps.empty()) { - cSteps.push_back(ChannelStep({0, 0}, end, start)); - std::sort(cSteps.begin(), cSteps.end()); - } - - std::vector cSegments; - cSegments.reserve(cSteps.size()); - - Bin2D currentBin = {bstart[0], bstart[1]}; - BinDelta2D lastDelta = {0, 0}; - Acts::Vector2 lastIntersect = start; - double lastPath = 0.; - for (auto& cStep : cSteps) { - currentBin[0] += lastDelta[0]; - currentBin[1] += lastDelta[1]; - double path = cStep.path - lastPath; - cSegments.push_back( - ChannelSegment(currentBin, {lastIntersect, cStep.intersect}, path)); - lastPath = cStep.path; - lastDelta = cStep.delta; - lastIntersect = cStep.intersect; - } - - return cSegments; -} diff --git a/Fatras/src/Digitization/Segmentizer.cpp b/Fatras/src/Digitization/Segmentizer.cpp new file mode 100644 index 00000000000..f1a70a9ab32 --- /dev/null +++ b/Fatras/src/Digitization/Segmentizer.cpp @@ -0,0 +1,210 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2020 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "ActsFatras/Digitization/Segmentizer.hpp" + +#include "Acts/Surfaces/Surface.hpp" +#include "Acts/Surfaces/detail/IntersectionHelper2D.hpp" +#include "Acts/Utilities/BinUtility.hpp" +#include "Acts/Utilities/BinningType.hpp" +#include "Acts/Utilities/Helpers.hpp" +#include "Acts/Utilities/Intersection.hpp" + +#include +#include +#include + +namespace { + +using ChannelStep = ActsFatras::Segmentizer::ChannelStep; +using ChannelSegment = ActsFatras::Segmentizer::ChannelSegment; +using Bin2D = ActsFatras::Segmentizer::Bin2D; +using BinDelta2D = ActsFatras::Segmentizer::BinDelta2D; + +auto stepsToSegments(std::vector& cSteps, + const Acts::Vector2& start, const Bin2D& bstart) { + std::vector cSegments; + cSegments.reserve(cSteps.size()); + + Bin2D currentBin = {bstart[0], bstart[1]}; + BinDelta2D lastDelta = {0, 0}; + Acts::Vector2 lastIntersect = start; + double lastPath = 0.; + + for (auto& cStep : cSteps) { + currentBin[0] += lastDelta[0]; + currentBin[1] += lastDelta[1]; + double path = cStep.path - lastPath; + cSegments.push_back( + ChannelSegment(currentBin, {lastIntersect, cStep.intersect}, path)); + lastPath = cStep.path; + lastDelta = cStep.delta; + lastIntersect = cStep.intersect; + } + + return cSegments; +} + +std::vector segmentPlaneSurface( + const Acts::BinUtility& segmentation, const Acts::Vector2& start, + const Acts::Vector2& end) { + // Full path length - the full channel + auto segment2d = (end - start); + Bin2D bstart = {0, 0}; + Bin2D bend = {0, 0}; + + // Get the segmentation and convert it to lines & arcs + bstart = {static_cast(segmentation.bin(start, 0)), + static_cast(segmentation.bin(start, 1))}; + bend = {static_cast(segmentation.bin(end, 0)), + static_cast(segmentation.bin(end, 1))}; + // Fast single channel exit + if (bstart == bend) { + return {ChannelSegment(bstart, {start, end}, segment2d.norm())}; + } + + std::vector cSteps; + + // The lines channel segment lines along x + if (bstart[0] != bend[0]) { + double k = segment2d.y() / segment2d.x(); + double d = start.y() - k * start.x(); + + const auto& xboundaries = segmentation.binningData()[0].boundaries(); + std::vector xbbounds = { + xboundaries.begin() + std::min(bstart[0], bend[0]) + 1, + xboundaries.begin() + std::max(bstart[0], bend[0]) + 1}; + for (const auto x : xbbounds) { + cSteps.push_back(ChannelStep{ + {(bstart[0] < bend[0] ? 1 : -1), 0}, {x, k * x + d}, start}); + } + } + // The lines channel segment lines along y + if (bstart[1] != bend[1]) { + double k = segment2d.x() / segment2d.y(); + double d = start.x() - k * start.y(); + const auto& yboundaries = segmentation.binningData()[1].boundaries(); + std::vector ybbounds = { + yboundaries.begin() + std::min(bstart[1], bend[1]) + 1, + yboundaries.begin() + std::max(bstart[1], bend[1]) + 1}; + for (const auto y : ybbounds) { + cSteps.push_back(ChannelStep{ + {0, (bstart[1] < bend[1] ? 1 : -1)}, {k * y + d, y}, start}); + } + } + + // Register the last step if successful + if (!cSteps.empty()) { + cSteps.push_back(ChannelStep({0, 0}, end, start)); + std::sort(cSteps.begin(), cSteps.end()); + } + + return stepsToSegments(cSteps, start, bstart); +} + +std::vector segmentDiscSurface( + const Acts::GeometryContext& geoCtx, const Acts::Surface& surface, + const Acts::BinUtility& segmentation, const Acts::Vector2& start, + const Acts::Vector2& end) { + // Full path length - the full channel + auto segment2d = (end - start); + Bin2D bstart = {0, 0}; + Bin2D bend = {0, 0}; + + Acts::Vector2 pstart(Acts::VectorHelpers::perp(start), + Acts::VectorHelpers::phi(start)); + Acts::Vector2 pend(Acts::VectorHelpers::perp(end), + Acts::VectorHelpers::phi(end)); + + // Get the segmentation and convert it to lines & arcs + bstart = {static_cast(segmentation.bin(pstart, 0)), + static_cast(segmentation.bin(pstart, 1))}; + bend = {static_cast(segmentation.bin(pend, 0)), + static_cast(segmentation.bin(pend, 1))}; + + // Fast single channel exit + if (bstart == bend) { + return {ChannelSegment(bstart, {start, end}, segment2d.norm())}; + } + + std::vector cSteps; + + double phistart = pstart[1]; + double phiend = pend[1]; + + // The radial boundaries + if (bstart[0] != bend[0]) { + const auto& rboundaries = segmentation.binningData()[0].boundaries(); + std::vector rbbounds = { + rboundaries.begin() + std::min(bstart[0], bend[0]) + 1, + rboundaries.begin() + std::max(bstart[0], bend[0]) + 1}; + for (const auto& r : rbbounds) { + auto radIntersection = + Acts::detail::IntersectionHelper2D::intersectCircleSegment( + r, std::min(phistart, phiend), std::max(phistart, phiend), start, + (end - start).normalized()); + cSteps.push_back(ChannelStep{{(bstart[0] < bend[0] ? 1 : -1), 0}, + radIntersection.position(), + start}); + } + } + // The phi boundaries + if (bstart[1] != bend[1]) { + double referenceR = surface.binningPositionValue(geoCtx, Acts::binR); + Acts::Vector2 origin = {0., 0.}; + const auto& phiboundaries = segmentation.binningData()[1].boundaries(); + std::vector phibbounds = { + phiboundaries.begin() + std::min(bstart[1], bend[1]) + 1, + phiboundaries.begin() + std::max(bstart[1], bend[1]) + 1}; + + for (const auto& phi : phibbounds) { + Acts::Vector2 philine(referenceR * std::cos(phi), + referenceR * std::sin(phi)); + auto phiIntersection = + Acts::detail::IntersectionHelper2D::intersectSegment( + origin, philine, start, (end - start).normalized()); + cSteps.push_back(ChannelStep{{0, (bstart[1] < bend[1] ? 1 : -1)}, + phiIntersection.position(), + start}); + } + } + + // Register the last step if successful + if (!cSteps.empty()) { + cSteps.push_back(ChannelStep({0, 0}, end, start)); + std::sort(cSteps.begin(), cSteps.end()); + } + + return stepsToSegments(cSteps, start, bstart); +} +} // namespace + +std::vector +ActsFatras::Segmentizer::segments(const Acts::GeometryContext& geoCtx, + const Acts::Surface& surface, + const Acts::BinUtility& segmentation, + const Segment2D& segment) const { + // Return if the segmentation is not two-dimensional + // (strips need to have one bin along the strip) + if (segmentation.dimensions() != 2) { + return {}; + } + + // Start and end point + const auto& start = segment[0]; + const auto& end = segment[1]; + + if (surface.type() == Acts::Surface::SurfaceType::Plane) { + return segmentPlaneSurface(segmentation, start, end); + } else if (surface.type() == Acts::Surface::SurfaceType::Disc) { + return segmentDiscSurface(geoCtx, surface, segmentation, start, end); + } else { + // TODO propagate error code + throw std::runtime_error("Cannot segment surface"); + } +} diff --git a/Plugins/ExaTrkX/CMakeLists.txt b/Plugins/ExaTrkX/CMakeLists.txt index 7202deeb508..ef180fd648a 100644 --- a/Plugins/ExaTrkX/CMakeLists.txt +++ b/Plugins/ExaTrkX/CMakeLists.txt @@ -1,6 +1,7 @@ set(SOURCES src/buildEdges.cpp src/ExaTrkXPipeline.cpp + src/CudaInfo.cpp ) if(ACTS_EXATRKX_ENABLE_ONNX) @@ -17,6 +18,7 @@ if(ACTS_EXATRKX_ENABLE_TORCH) src/TorchMetricLearning.cpp src/BoostTrackBuilding.cpp src/TorchTruthGraphMetricsHook.cpp + src/TorchGraphStoreHook.cpp ) endif() @@ -102,6 +104,15 @@ if(ACTS_EXATRKX_ENABLE_TORCH) endif() +# add_executable(EdgeClassifierShiftChecker src/EdgeClassifierShiftChecker.cpp) +# target_link_libraries(EdgeClassifierShiftChecker PUBLIC ActsPluginExaTrkX ${TORCH_LIBRARIES}) +# +# add_executable(GraphConstructorShiftChecker src/GraphConstructorShiftChecker.cpp) +# target_link_libraries(GraphConstructorShiftChecker PUBLIC ActsPluginExaTrkX ${TORCH_LIBRARIES}) +# +# add_executable(CheckPipelineShift.cpp src/CheckPipelineShift.cpp) +# target_link_libraries(CheckPipelineShift.cpp PUBLIC ActsPluginExaTrkX ${TORCH_LIBRARIES}) + install( TARGETS ActsPluginExaTrkX EXPORT ActsPluginExaTrkXTargets diff --git a/Plugins/ExaTrkX/include/Acts/Plugins/ExaTrkX/ExaTrkXPipeline.hpp b/Plugins/ExaTrkX/include/Acts/Plugins/ExaTrkX/ExaTrkXPipeline.hpp index e6810eb22ef..19d3a5db985 100644 --- a/Plugins/ExaTrkX/include/Acts/Plugins/ExaTrkX/ExaTrkXPipeline.hpp +++ b/Plugins/ExaTrkX/include/Acts/Plugins/ExaTrkX/ExaTrkXPipeline.hpp @@ -31,8 +31,9 @@ struct ExaTrkXTiming { class ExaTrkXHook { public: - virtual ~ExaTrkXHook() {} - virtual void operator()(const std::any &, const std::any &) const {}; + virtual ~ExaTrkXHook(){}; + virtual void operator()(const std::any &, const std::any &, + const std::any &) const {}; }; class ExaTrkXPipeline { diff --git a/Plugins/ExaTrkX/include/Acts/Plugins/ExaTrkX/Stages.hpp b/Plugins/ExaTrkX/include/Acts/Plugins/ExaTrkX/Stages.hpp index 048f56bfd3c..a6107d9c67b 100644 --- a/Plugins/ExaTrkX/include/Acts/Plugins/ExaTrkX/Stages.hpp +++ b/Plugins/ExaTrkX/include/Acts/Plugins/ExaTrkX/Stages.hpp @@ -11,6 +11,8 @@ #include #include +#include + namespace Acts { // TODO maybe replace std::any with some kind of variant, diff --git a/Plugins/ExaTrkX/include/Acts/Plugins/ExaTrkX/TorchGraphStoreHook.hpp b/Plugins/ExaTrkX/include/Acts/Plugins/ExaTrkX/TorchGraphStoreHook.hpp new file mode 100644 index 00000000000..172b96fce20 --- /dev/null +++ b/Plugins/ExaTrkX/include/Acts/Plugins/ExaTrkX/TorchGraphStoreHook.hpp @@ -0,0 +1,34 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2023 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/Plugins/ExaTrkX/ExaTrkXPipeline.hpp" +#include "Acts/Plugins/ExaTrkX/detail/CantorEdge.hpp" +#include "Acts/Utilities/Logger.hpp" + +namespace Acts { + +class TorchGraphStoreHook : public ExaTrkXHook { + public: + using Graph = std::pair, std::vector>; + + private: + std::unique_ptr m_storedGraph; + + public: + TorchGraphStoreHook(); + ~TorchGraphStoreHook() override {} + + void operator()(const std::any &, const std::any &edges, + const std::any &weights) const override; + + const Graph &storedGraph() const { return *m_storedGraph; } +}; + +} // namespace Acts diff --git a/Plugins/ExaTrkX/include/Acts/Plugins/ExaTrkX/TorchTruthGraphMetricsHook.hpp b/Plugins/ExaTrkX/include/Acts/Plugins/ExaTrkX/TorchTruthGraphMetricsHook.hpp index a13c9de984d..f971ae2992e 100644 --- a/Plugins/ExaTrkX/include/Acts/Plugins/ExaTrkX/TorchTruthGraphMetricsHook.hpp +++ b/Plugins/ExaTrkX/include/Acts/Plugins/ExaTrkX/TorchTruthGraphMetricsHook.hpp @@ -25,7 +25,8 @@ class TorchTruthGraphMetricsHook : public ExaTrkXHook { std::unique_ptr l); ~TorchTruthGraphMetricsHook() override {} - void operator()(const std::any &, const std::any &edges) const override; + void operator()(const std::any &, const std::any &edges, + const std::any &) const override; }; } // namespace Acts diff --git a/Plugins/ExaTrkX/include/Acts/Plugins/ExaTrkX/detail/BoostTrackBuildingUtils.hpp b/Plugins/ExaTrkX/include/Acts/Plugins/ExaTrkX/detail/BoostTrackBuildingUtils.hpp new file mode 100644 index 00000000000..143a9ae489b --- /dev/null +++ b/Plugins/ExaTrkX/include/Acts/Plugins/ExaTrkX/detail/BoostTrackBuildingUtils.hpp @@ -0,0 +1,110 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2023 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/Utilities/Logger.hpp" + +#include +#include +#include +#include + +namespace Acts::detail { + +struct SubgraphNodePredicate { + std::size_t subgraphId = 0; + std::vector *labels = nullptr; + + template + bool operator()(const vertex_t &v) const { + return labels->at(v) == subgraphId; + } +}; + +template +struct SubgraphEdgePredicate { + graph_t *g = nullptr; + template + bool operator()(const edge_t &e) const { + return (*g)[e].weight > 0.0; + } +}; + +template +bool isCleanSubgraph(const graph_t &graph) { + std::size_t nStart = 0; + std::size_t nStop = 0; + + for (const auto &n : boost::make_iterator_range(boost::vertices(graph))) { + nStart += (boost::in_degree(n, graph) == 0); + nStop += (boost::out_degree(n, graph) == 0); + } + + return nStart == 1 and nStop == 1; +} + +template +void cleanSubgraphs(graph_t &graph, + const Acts::Logger &logger = Acts::getDummyLogger()) { + SubgraphEdgePredicate edgeFilter{&graph}; + + using Subgraph = + boost::filtered_graph, + SubgraphNodePredicate>; + + std::vector connectedComponentLabels(boost::num_vertices(graph)); + auto nSubgraphs = + boost::connected_components(graph, connectedComponentLabels.data()); + + for (auto i = 0ul; i < nSubgraphs; ++i) { + SubgraphNodePredicate nodeFilter{i, &connectedComponentLabels}; + Subgraph subgraph(graph, edgeFilter, nodeFilter); + + if (isCleanSubgraph(subgraph)) { + continue; + } + + while (true) { + // Find edge with minium weight (edges with weight 0 should not occur + // because of edge filter) However, we only check edges which are + // branching edges + auto [edgeBegin, edgeEnd] = boost::edges(subgraph); + float minWeight = std::numeric_limits::max(); + typename Subgraph::edge_iterator minEdgeDesc = edgeEnd; + for (auto it = edgeBegin; it != edgeEnd; ++it) { + assert(graph[*it].weight > 0); + if (boost::out_degree(boost::source(*it, subgraph), subgraph) < 2 && + boost::in_degree(boost::target(*it, subgraph), subgraph) < 2) { + continue; + } + if (graph[*it].weight < minWeight) { + minWeight = graph[*it].weight; + minEdgeDesc = it; + } + } + + // All edges are clean + if (minEdgeDesc == edgeEnd) { + break; + } + + // Set edge to 0 (should effectively remove the edge from the filtered + // graph) + ACTS_VERBOSE("Remove edge " << boost::source(*minEdgeDesc, graph) << ", " + << boost::target(*minEdgeDesc, graph)); + graph[*minEdgeDesc].weight = 0.0; + } + } + + // Finally remove all edges with weight == 0 + boost::remove_edge_if([&](auto ed) { return graph[ed].weight == 0.f; }, + graph); +} + +} // namespace Acts::detail diff --git a/Plugins/ExaTrkX/include/Acts/Plugins/ExaTrkX/detail/CudaInfo.hpp b/Plugins/ExaTrkX/include/Acts/Plugins/ExaTrkX/detail/CudaInfo.hpp new file mode 100644 index 00000000000..f50081b9eaa --- /dev/null +++ b/Plugins/ExaTrkX/include/Acts/Plugins/ExaTrkX/detail/CudaInfo.hpp @@ -0,0 +1,17 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2022 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +#include +#include + +namespace Acts::detail { +/// Returns device count +std::size_t cudaNumDevices(); +} // namespace Acts::detail diff --git a/Plugins/ExaTrkX/src/BoostTrackBuilding.cpp b/Plugins/ExaTrkX/src/BoostTrackBuilding.cpp index 01e80930960..f49b0063ad4 100644 --- a/Plugins/ExaTrkX/src/BoostTrackBuilding.cpp +++ b/Plugins/ExaTrkX/src/BoostTrackBuilding.cpp @@ -15,33 +15,10 @@ #include #include #include +#include #include -namespace { -template -auto weaklyConnectedComponents(vertex_t numNodes, - boost::beast::span& rowIndices, - boost::beast::span& colIndices, - boost::beast::span& edgeWeights, - std::vector& trackLabels) { - using Graph = - boost::adjacency_list; - - Graph g(numNodes); - - for (const auto [row, col, weight] : - Acts::zip(rowIndices, colIndices, edgeWeights)) { - boost::add_edge(row, col, weight, g); - } - - return boost::connected_components(g, &trackLabels[0]); -} -} // namespace +using namespace torch::indexing; namespace Acts { @@ -49,7 +26,15 @@ std::vector> BoostTrackBuilding::operator()( std::any nodes, std::any edges, std::any weights, std::vector& spacepointIDs, int) { ACTS_DEBUG("Start track building"); + + // Get nodes + const auto nodeTensor = std::any_cast(nodes).to(torch::kCPU); + assert(static_cast(nodeTensor.size(0)) == spacepointIDs.size()); + + // Get edges const auto edgeTensor = std::any_cast(edges).to(torch::kCPU); + + // Get weights const auto edgeWeightTensor = std::any_cast(weights).to(torch::kCPU); @@ -74,11 +59,42 @@ std::vector> BoostTrackBuilding::operator()( boost::beast::span edgeWeights(edgeWeightTensor.data_ptr(), numEdges); - std::vector trackLabels(numSpacepoints); + // TODO is this clone necessary? + const auto radiusTensor = nodeTensor.index({Slice{}, 0}).clone(); + boost::beast::span nodeRadius(radiusTensor.data_ptr(), + radiusTensor.numel()); + + // Construct Graph + struct EdgeProperty { + weight_t weight; + }; - auto numberLabels = weaklyConnectedComponents( - numSpacepoints, rowIndices, colIndices, edgeWeights, trackLabels); + using Graph = + boost::adjacency_list; + + Graph g(numSpacepoints); + + for (const auto [row, col, weight] : + Acts::zip(rowIndices, colIndices, edgeWeights)) { + const auto rowRadius = *(nodeRadius.begin() + row); + const auto colRadius = *(nodeRadius.begin() + col); + if (rowRadius < colRadius) { + boost::add_edge(row, col, EdgeProperty{weight}, g); + } else { + boost::add_edge(col, row, EdgeProperty{weight}, g); + } + } + + // Make final connected components + std::vector trackLabels(numSpacepoints); + auto numberLabels = boost::connected_components(g, &trackLabels[0]); + // Label edges ACTS_VERBOSE("Number of track labels: " << trackLabels.size()); ACTS_VERBOSE("Number of unique track labels: " << [&]() { std::vector sorted(trackLabels); diff --git a/Plugins/ExaTrkX/src/CheckPipelineShift.cpp b/Plugins/ExaTrkX/src/CheckPipelineShift.cpp new file mode 100644 index 00000000000..c8c2bff19d5 --- /dev/null +++ b/Plugins/ExaTrkX/src/CheckPipelineShift.cpp @@ -0,0 +1,107 @@ + +#include +#include +#include +#include + +#include +#include + +#include +#include + +const auto device = torch::cuda::is_available() ? torch::kCUDA : torch::kCPU; + +auto run_pipeline(std::vector &data, std::vector &spacepointIDs, + const std::string &path, + Acts::Logging::Level l = Acts::Logging::VERBOSE) { + Acts::TorchMetricLearning::Config gcCfg; + gcCfg.knnVal = 100; + gcCfg.rVal = 0.1; + gcCfg.modelPath = path + "/embedding.pt"; + gcCfg.numFeatures = 7; + + Acts::TorchEdgeClassifier::Config fltCfg; + fltCfg.modelPath = path + "/filter.pt"; + fltCfg.nChunks = 0; + fltCfg.cut = 0.5; + fltCfg.numFeatures = 3; + fltCfg.undirected = false; + + Acts::TorchEdgeClassifier::Config gnnCfg; + gnnCfg.modelPath = path + "/gnn.pt"; + gnnCfg.nChunks = 0; + gnnCfg.cut = 0.01; + gnnCfg.numFeatures = 3; + gnnCfg.undirected = true; + + auto gc = std::make_shared( + gcCfg, Acts::getDefaultLogger("gc", l)); + auto flt = std::make_shared( + fltCfg, Acts::getDefaultLogger("flt", l)); + auto gnn = std::make_shared( + gnnCfg, Acts::getDefaultLogger("gnn", l)); + auto trk = std::make_shared( + Acts::getDefaultLogger("trk", l)); + + Acts::Pipeline pipeline(gc, {flt, gnn}, trk, + Acts::getDefaultLogger("pipeline", l)); + + return pipeline.run(data, spacepointIDs); +} + +void checkShiftInvariance(const std::string &path) { + auto features = torch::rand({100, 7}).to(torch::kFloat); + std::vector feature_vec(features.data_ptr(), + features.data_ptr() + features.numel()); + + std::vector spacepointIDs(100); + std::iota(spacepointIDs.begin(), spacepointIDs.end(), 0); + + // Reference + auto tracks1 = run_pipeline(feature_vec, spacepointIDs, path); + + // Shift + std::rotate(spacepointIDs.begin(), spacepointIDs.begin() + 1, + spacepointIDs.end()); + auto features_rolled = torch::roll(features, 1, 0).clone(); + std::vector feature_rolled_vec( + features_rolled.data_ptr(), + features_rolled.data_ptr() + features_rolled.numel()); + + // Shifted run + auto tracks2 = run_pipeline(feature_rolled_vec, spacepointIDs, path); + + // Print + std::vector sizes1(tracks1.size()), sizes2(tracks2.size()); + std::transform(tracks1.begin(), tracks1.end(), sizes1.begin(), + [](const auto &t) { return t.size(); }); + std::transform(tracks2.begin(), tracks2.end(), sizes2.begin(), + [](const auto &t) { return t.size(); }); + + std::sort(sizes1.begin(), sizes1.end()); + std::sort(sizes2.begin(), sizes2.end()); + + std::cout << "Tracks 1:"; + std::copy(sizes1.begin(), sizes1.end(), + std::ostream_iterator(std::cout, " ")); + std::cout << std::endl; + + std::cout << "Tracks 2:"; + std::copy(sizes2.begin(), sizes2.end(), + std::ostream_iterator(std::cout, " ")); + std::cout << std::endl; +} + +int main(int argc, char **argv) { + std::vector args(argv, argv + argc); + + std::cout << "check model only:\n"; + checkShiftInvariance(args.at(1)); + + // std::cout << "check stage:\n"; + // checkClassifierStage(args.at(1)); + // + // std::cout << "check pipeline:\n"; + // checkPipeline(args.at(1)); +} diff --git a/Plugins/ExaTrkX/src/CudaInfo.cpp b/Plugins/ExaTrkX/src/CudaInfo.cpp new file mode 100644 index 00000000000..fdce8cbe28c --- /dev/null +++ b/Plugins/ExaTrkX/src/CudaInfo.cpp @@ -0,0 +1,11 @@ +#include "Acts/Plugins/ExaTrkX/detail/CudaInfo.hpp" + +#include "torch/torch.h" + +std::size_t Acts::detail::cudaNumDevices() { + if (not torch::cuda::is_available()) { + return 0; + } + + return torch::cuda::device_count(); +} diff --git a/Plugins/ExaTrkX/src/EdgeClassifierShiftChecker.cpp b/Plugins/ExaTrkX/src/EdgeClassifierShiftChecker.cpp new file mode 100644 index 00000000000..e6e40ae2287 --- /dev/null +++ b/Plugins/ExaTrkX/src/EdgeClassifierShiftChecker.cpp @@ -0,0 +1,149 @@ + +#include +#include + +#include +#include + +#include +#include + +const auto device = torch::cuda::is_available() ? torch::kCUDA : torch::kCPU; + +struct DummyGraphConstruction : Acts::GraphConstructionBase { + torch::Tensor features; + torch::Tensor edges; + + std::tuple operator()(std::vector &, std::size_t, + int) override { + return {features, edges}; + } +}; + +struct DummyTrackBuilder : Acts::TrackBuildingBase { + torch::Tensor out_features; + torch::Tensor out_edges; + torch::Tensor out_weights; + + std::vector> operator()(std::any nodes, std::any edges, + std::any edgeWeights, + std::vector &, int) override { + out_features = std::any_cast(nodes); + out_edges = std::any_cast(edges); + out_weights = std::any_cast(edgeWeights); + + return {}; + } +}; + +std::tuple getTensors(int nNodes = 100, + int nEdges = 20) { + return {torch::rand({nNodes, 3}).to(torch::kFloat32).to(device), + torch::randint(0, nNodes, {2, nEdges}).to(device)}; +} + +std::tuple modifyTensors(const torch::Tensor &features, + const torch::Tensor &edges) { + return {torch::roll(features, 1, 0).clone(), + ((edges + 1) % features.size(0)).clone()}; +} + +void printDiff(const at::Tensor &out1, const at::Tensor &out2) { + auto diff = torch::abs(out1 - out2).to(torch::kCPU); + + std::cout << "diff: "; + std::copy(diff.data_ptr(), diff.data_ptr() + diff.numel(), + std::ostream_iterator(std::cout, " ")); + std::cout << std::endl; +} + +void checkModelShiftInvariance(const std::string &path) { + auto model = torch::jit::load(path); + model.to(device); + + auto [features, edges] = getTensors(); + + std::vector input; + input.push_back(features); + input.push_back(edges); + + auto output = model.forward(input).toTensor(); + output = std::get<0>(torch::sort(output)); + + auto [features2, edges2] = modifyTensors(features, edges); + + input.clear(); + input.push_back(features2); + input.push_back(edges2); + + auto output2 = model.forward(input).toTensor(); + output2 = std::get<0>(torch::sort(output2)); + + printDiff(output, output2); +} + +void checkClassifierStage(const std::string &path) { + Acts::TorchEdgeClassifier::Config cfg; + cfg.modelPath = path; + cfg.nChunks = 1; + cfg.cut = 0.0; + cfg.numFeatures = 3; + cfg.undirected = false; + + auto logger = Acts::getDefaultLogger("test", Acts::Logging::INFO); + Acts::TorchEdgeClassifier clf(cfg, std::move(logger)); + + auto [features, edges] = getTensors(); + auto [nodes_out, features_out, output] = clf(features, edges); + + auto [features2, edges2] = modifyTensors(features, edges); + auto [nodes_out2, features_out2, output2] = clf(features, edges); + + auto output_tensor = std::any_cast(output); + auto output_tensor2 = std::any_cast(output2); + + printDiff(output_tensor, output_tensor2); +} + +void checkPipeline(const std::string &path) { + Acts::TorchEdgeClassifier::Config cfg; + cfg.modelPath = path; + cfg.nChunks = 1; + cfg.cut = 0.0; + cfg.numFeatures = 3; + cfg.undirected = false; + + auto gc = std::make_shared(); + auto cls = std::make_shared( + cfg, Acts::getDefaultLogger("test", Acts::Logging::INFO)); + auto trk = std::make_shared(); + + Acts::Pipeline pipeline(gc, {cls}, trk, + Acts::getDefaultLogger("test", Acts::Logging::INFO)); + + std::vector dummyData; + std::vector dummyIds; + + std::tie(gc->features, gc->edges) = getTensors(); + pipeline.run(dummyData, dummyIds); + auto output = std::get<0>(torch::sort(trk->out_weights.clone())); + + std::tie(gc->features, gc->edges) = modifyTensors(gc->features, gc->edges); + pipeline.run(dummyData, dummyIds); + auto output2 = std::get<0>(torch::sort(trk->out_weights.clone())); + + printDiff(output, output2); +} + +int main(int argc, char **argv) { + std::vector args(argv, argv + argc); + + std::cout << "check model only:\n"; + checkModelShiftInvariance(args.at(1)); + + std::cout << "check stage:\n"; + checkClassifierStage(args.at(1)); + + std::cout << "check pipeline:\n"; + checkPipeline(args.at(1)); +} diff --git a/Plugins/ExaTrkX/src/ExaTrkXPipeline.cpp b/Plugins/ExaTrkX/src/ExaTrkXPipeline.cpp index 8c408413c16..3f8e88150f0 100644 --- a/Plugins/ExaTrkX/src/ExaTrkXPipeline.cpp +++ b/Plugins/ExaTrkX/src/ExaTrkXPipeline.cpp @@ -44,7 +44,7 @@ std::vector> ExaTrkXPipeline::run( timing->graphBuildingTime = t1 - t0; } - hook(nodes, edges); + hook(nodes, edges, {}); std::any edge_weights; timing->classifierTimes.clear(); @@ -63,7 +63,7 @@ std::vector> ExaTrkXPipeline::run( edges = std::move(newEdges); edge_weights = std::move(newWeights); - hook(nodes, edges); + hook(nodes, edges, edge_weights); } t0 = std::chrono::high_resolution_clock::now(); diff --git a/Plugins/ExaTrkX/src/GraphConstructorShiftChecker.cpp b/Plugins/ExaTrkX/src/GraphConstructorShiftChecker.cpp new file mode 100644 index 00000000000..26187bf0eb7 --- /dev/null +++ b/Plugins/ExaTrkX/src/GraphConstructorShiftChecker.cpp @@ -0,0 +1,73 @@ + +#include +#include + +#include +#include + +#include +#include + +const auto device = torch::cuda::is_available() ? torch::kCUDA : torch::kCPU; + +void printDiff(const at::Tensor &out1, const at::Tensor &out2) { + auto diff = torch::abs(out1 - out2).to(torch::kCPU); + + std::cout << "diff: "; + std::copy(diff.data_ptr(), diff.data_ptr() + diff.numel(), + std::ostream_iterator(std::cout, " ")); + std::cout << std::endl; +} + +void checkModelShiftInvariance(const std::string &path) { + auto model = torch::jit::load(path); + model.to(device); + + Acts::TorchMetricLearning::Config cfg; + cfg.knnVal = 100; + cfg.rVal = 0.1; + cfg.modelPath = path; + cfg.numFeatures = 7; + + Acts::TorchMetricLearning gc( + cfg, Acts::getDefaultLogger("test", Acts::Logging::INFO)); + + // Reference run + auto features = torch::rand({100, 7}).to(torch::kFloat); + // std::cout << features << std::endl; + std::vector feature_vec(features.data_ptr(), + features.data_ptr() + features.numel()); + + auto [nodes, edges] = gc(feature_vec, 100); + auto edges1Tensor = std::any_cast(edges); + edges1Tensor = std::get<0>(torch::sort(edges1Tensor, 0)); + + // Shifted run + auto features_rolled = torch::roll(features, 1, 0).clone(); + // std::cout << features_rolled << std::endl; + std::vector feature_rolled_vec( + features_rolled.data_ptr(), + features_rolled.data_ptr() + features_rolled.numel()); + + auto [nodes2, edges2] = gc(feature_rolled_vec, 100); + auto edges2Tensor = std::any_cast(edges2); + edges2Tensor = std::get<0>(torch::sort(edges2Tensor, 0)); + + // Print + auto shift = (edges1Tensor + 1) % 100; + std::cout << shift << std::endl; + std::cout << edges2Tensor << std::endl; +} + +int main(int argc, char **argv) { + std::vector args(argv, argv + argc); + + std::cout << "check model only:\n"; + checkModelShiftInvariance(args.at(1)); + + // std::cout << "check stage:\n"; + // checkClassifierStage(args.at(1)); + // + // std::cout << "check pipeline:\n"; + // checkPipeline(args.at(1)); +} diff --git a/Plugins/ExaTrkX/src/Pipeline.cpp b/Plugins/ExaTrkX/src/Pipeline.cpp new file mode 100644 index 00000000000..dd36cf7bb57 --- /dev/null +++ b/Plugins/ExaTrkX/src/Pipeline.cpp @@ -0,0 +1,60 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2023 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "Acts/Plugins/ExaTrkX/Pipeline.hpp" + +namespace Acts { + +Pipeline::Pipeline( + std::shared_ptr graphConstructor, + std::vector> edgeClassifiers, + std::shared_ptr trackBuilder, + std::unique_ptr logger) + : m_logger(std::move(logger)), + m_graphConstructor(graphConstructor), + m_edgeClassifiers(edgeClassifiers), + m_trackBuilder(trackBuilder) { + if (!m_graphConstructor) { + throw std::invalid_argument("Missing graph construction module"); + } + if (!m_trackBuilder) { + throw std::invalid_argument("Missing track building module"); + } + if (m_edgeClassifiers.empty() or + not std::all_of(m_edgeClassifiers.begin(), m_edgeClassifiers.end(), + [](const auto &a) { return static_cast(a); })) { + throw std::invalid_argument("Missing graph construction module"); + } +} + +std::vector> Pipeline::run(std::vector &features, + std::vector &spacepointIDs, + int deviceHint, + const PipelineHook &hook) const { + auto [nodes, edges] = + (*m_graphConstructor)(features, spacepointIDs.size(), deviceHint); + + hook(nodes, edges); + + std::any edge_weights; + + for (auto edgeClassifier : m_edgeClassifiers) { + auto [newNodes, newEdges, newWeights] = + (*edgeClassifier)(std::move(nodes), std::move(edges), deviceHint); + nodes = std::move(newNodes); + edges = std::move(newEdges); + edge_weights = std::move(newWeights); + + hook(nodes, edges); + } + + return (*m_trackBuilder)(std::move(nodes), std::move(edges), + std::move(edge_weights), spacepointIDs, deviceHint); +} + +} // namespace Acts diff --git a/Plugins/ExaTrkX/src/TorchEdgeClassifier.cpp b/Plugins/ExaTrkX/src/TorchEdgeClassifier.cpp index 8e9a1c46d7b..9a276869112 100644 --- a/Plugins/ExaTrkX/src/TorchEdgeClassifier.cpp +++ b/Plugins/ExaTrkX/src/TorchEdgeClassifier.cpp @@ -51,40 +51,73 @@ std::tuple TorchEdgeClassifier::operator()( auto nodes = std::any_cast(inputNodes).to(device); auto edgeList = std::any_cast(inputEdges).to(device); + auto model = m_model->clone(); + model.to(device); + if (m_cfg.numFeatures > nodes.size(1)) { throw std::runtime_error("requested more features then available"); } - std::vector results; - results.reserve(m_cfg.nChunks); - - auto edgeListTmp = - m_cfg.undirected ? torch::cat({edgeList, edgeList.flip(0)}, 1) : edgeList; - - std::vector inputTensors(2); - inputTensors[0] = m_cfg.numFeatures < nodes.size(1) - ? nodes.index({Slice{}, Slice{None, m_cfg.numFeatures}}) - : nodes; - - const auto chunks = at::chunk(at::arange(edgeListTmp.size(1)), m_cfg.nChunks); - for (const auto& chunk : chunks) { - ACTS_VERBOSE("Process chunk"); - inputTensors[1] = edgeListTmp.index({Slice(), chunk}); - - results.push_back(m_model->forward(inputTensors).toTensor()); - results.back().squeeze_(); - results.back().sigmoid_(); + torch::Tensor output; + + // Scope this to keep inference objects seperate + { + auto edgeListTmp = m_cfg.undirected + ? torch::cat({edgeList, edgeList.flip(0)}, 1) + : edgeList; + + std::vector inputTensors(2); + inputTensors[0] = + m_cfg.numFeatures < nodes.size(1) + ? nodes.index({Slice{}, Slice{None, m_cfg.numFeatures}}) + : nodes; + + if (m_cfg.nChunks > 1) { + std::vector results; + results.reserve(m_cfg.nChunks); + + auto chunks = at::chunk(edgeListTmp, m_cfg.nChunks, 1); + for (auto& chunk : chunks) { + ACTS_VERBOSE("Process chunk with shape" << chunk.sizes()); + inputTensors[1] = chunk; + + results.push_back(model.forward(inputTensors).toTensor()); + results.back().squeeze_(); + } + + output = torch::cat(results); + } else { + inputTensors[1] = edgeListTmp; + output = model.forward(inputTensors).toTensor(); + output.squeeze_(); + } } - auto output = torch::cat(results); + output.sigmoid_(); if (m_cfg.undirected) { - output = output.index({Slice(None, output.size(0) / 2)}); + auto newSize = output.size(0) / 2; + output = output.index({Slice(None, newSize)}); } ACTS_VERBOSE("Size after classifier: " << output.size(0)); +#if 0 + ACTS_VERBOSE("Slice of classified output:" << [&]() { + std::stringstream ss; + auto idxs = torch::argsort(output).to(torch::kInt64); + for (int i : {0, 1, static_cast(idxs.numel() / 2), -2, -1}) { + auto ii = idxs[i].item(); + ss << "\n" + << edgeList[0][ii].item() << ", " + << edgeList[1][ii].item() << " -> " + << output[ii].item(); + } + return ss.str(); + }()); +#else ACTS_VERBOSE("Slice of classified output:\n" << output.slice(/*dim=*/0, /*start=*/0, /*end=*/9)); +#endif printCudaMemInfo(logger()); torch::Tensor mask = output > m_cfg.cut; diff --git a/Plugins/ExaTrkX/src/TorchGraphStoreHook.cpp b/Plugins/ExaTrkX/src/TorchGraphStoreHook.cpp new file mode 100644 index 00000000000..0d5bd6b8c65 --- /dev/null +++ b/Plugins/ExaTrkX/src/TorchGraphStoreHook.cpp @@ -0,0 +1,33 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2023 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include "Acts/Plugins/ExaTrkX/TorchGraphStoreHook.hpp" + +#include "Acts/Plugins/ExaTrkX/detail/TensorVectorConversion.hpp" + +#include + +Acts::TorchGraphStoreHook::TorchGraphStoreHook() { + m_storedGraph = std::make_unique(); +} + +void Acts::TorchGraphStoreHook::operator()(const std::any&, + const std::any& edges, + const std::any& weights) const { + if (not weights.has_value()) { + return; + } + + m_storedGraph->first = detail::tensor2DToVector( + std::any_cast(edges).t()); + + auto cpuWeights = std::any_cast(weights).to(torch::kCPU); + m_storedGraph->second = + std::vector(cpuWeights.data_ptr(), + cpuWeights.data_ptr() + cpuWeights.numel()); +} diff --git a/Plugins/ExaTrkX/src/TorchTruthGraphMetricsHook.cpp b/Plugins/ExaTrkX/src/TorchTruthGraphMetricsHook.cpp index 50712797c39..851142274b9 100644 --- a/Plugins/ExaTrkX/src/TorchTruthGraphMetricsHook.cpp +++ b/Plugins/ExaTrkX/src/TorchTruthGraphMetricsHook.cpp @@ -46,7 +46,8 @@ Acts::TorchTruthGraphMetricsHook::TorchTruthGraphMetricsHook( } void Acts::TorchTruthGraphMetricsHook::operator()(const std::any&, - const std::any& edges) const { + const std::any& edges, + const std::any&) const { // We need to transpose the edges here for the right memory layout const auto edgeIndex = Acts::detail::tensor2DToVector( std::any_cast(edges).t()); diff --git a/Plugins/ExaTrkX/src/buildEdges.cpp b/Plugins/ExaTrkX/src/buildEdges.cpp index 1b133291d7f..b31047eca48 100644 --- a/Plugins/ExaTrkX/src/buildEdges.cpp +++ b/Plugins/ExaTrkX/src/buildEdges.cpp @@ -86,19 +86,19 @@ torch::Tensor Acts::detail::buildEdgesFRNN(torch::Tensor &embedFeatures, int G = -1; // Set up grid properties - torch::Tensor grid_min; - torch::Tensor grid_max; - torch::Tensor grid_size; + at::Tensor grid_min; + at::Tensor grid_max; + at::Tensor grid_size; - torch::Tensor embedTensor = embedFeatures.reshape({1, numSpacepoints, dim}); - torch::Tensor gridParamsCuda = + at::Tensor embedTensor = embedFeatures.reshape({1, numSpacepoints, dim}); + at::Tensor gridParamsCuda = torch::zeros({batch_size, grid_params_size}, device).to(torch::kFloat32); - torch::Tensor r_tensor = torch::full({batch_size}, rVal, device); - torch::Tensor lengths = torch::full({batch_size}, numSpacepoints, device); + at::Tensor r_tensor = torch::full({batch_size}, rVal, device); + at::Tensor lengths = torch::full({batch_size}, numSpacepoints, device); // build the grid for (int i = 0; i < batch_size; i++) { - torch::Tensor allPoints = + at::Tensor allPoints = embedTensor.index({i, Slice(None, lengths.index({i}).item().to()), Slice(None, grid_dim)}); grid_min = std::get<0>(allPoints.min(0)); @@ -129,28 +129,28 @@ torch::Tensor Acts::detail::buildEdgesFRNN(torch::Tensor &embedFeatures, } } - torch::Tensor pc_grid_cnt = + at::Tensor pc_grid_cnt = torch::zeros({batch_size, G}, device).to(torch::kInt32); - torch::Tensor pc_grid_cell = + at::Tensor pc_grid_cell = torch::full({batch_size, numSpacepoints}, -1, device).to(torch::kInt32); - torch::Tensor pc_grid_idx = + at::Tensor pc_grid_idx = torch::full({batch_size, numSpacepoints}, -1, device).to(torch::kInt32); // put spacepoints into the grid InsertPointsCUDA(embedTensor, lengths.to(torch::kInt64), gridParamsCuda, pc_grid_cnt, pc_grid_cell, pc_grid_idx, G); - torch::Tensor pc_grid_off = + at::Tensor pc_grid_off = torch::full({batch_size, G}, 0, device).to(torch::kInt32); - torch::Tensor grid_params = gridParamsCuda.to(torch::kCPU); + at::Tensor grid_params = gridParamsCuda.to(torch::kCPU); // for loop seems not to be necessary anymore pc_grid_off = PrefixSumCUDA(pc_grid_cnt, grid_params); - torch::Tensor sorted_points = + at::Tensor sorted_points = torch::zeros({batch_size, numSpacepoints, dim}, device) .to(torch::kFloat32); - torch::Tensor sorted_points_idxs = + at::Tensor sorted_points_idxs = torch::full({batch_size, numSpacepoints}, -1, device).to(torch::kInt32); CountingSortCUDA(embedTensor, lengths.to(torch::kInt64), pc_grid_cell, @@ -163,9 +163,9 @@ torch::Tensor Acts::detail::buildEdgesFRNN(torch::Tensor &embedFeatures, gridParamsCuda.to(torch::kFloat32), kVal, r_tensor, r_tensor * r_tensor); torch::Tensor positiveIndices = indices >= 0; - torch::Tensor repeatRange = torch::arange(positiveIndices.size(1), device) - .repeat({1, positiveIndices.size(2), 1}) - .transpose(1, 2); + at::Tensor repeatRange = torch::arange(positiveIndices.size(1), device) + .repeat({1, positiveIndices.size(2), 1}) + .transpose(1, 2); torch::Tensor stackedEdges = torch::stack( {repeatRange.index({positiveIndices}), indices.index({positiveIndices})}); diff --git a/Tests/UnitTests/Examples/Algorithms/Digitization/ModuleClustersTests.cpp b/Tests/UnitTests/Examples/Algorithms/Digitization/ModuleClustersTests.cpp index 5503d2d466c..9fb54404dbb 100644 --- a/Tests/UnitTests/Examples/Algorithms/Digitization/ModuleClustersTests.cpp +++ b/Tests/UnitTests/Examples/Algorithms/Digitization/ModuleClustersTests.cpp @@ -10,7 +10,7 @@ #include "Acts/Utilities/BinningData.hpp" #include "ActsExamples/Digitization/ModuleClusters.hpp" -#include "ActsFatras/Digitization/Channelizer.hpp" +#include "ActsFatras/Digitization/Segmentizer.hpp" using namespace Acts; using namespace ActsFatras; @@ -23,9 +23,9 @@ DigitizedParameters makeDigitizationParameters(const Vector2 &position, const BinUtility &binUtility) { auto [binX, binY, _] = binUtility.binTriple((Vector3() << position, 0).finished()); - Channelizer::Bin2D bin = {(Channelizer::Bin2D::value_type)binX, - (Channelizer::Bin2D::value_type)binY}; - Channelizer::Segment2D segment = {position, position}; + Segmentizer::Bin2D bin = {(Segmentizer::Bin2D::value_type)binX, + (Segmentizer::Bin2D::value_type)binY}; + Segmentizer::Segment2D segment = {position, position}; double activation = 1; Cluster::Cell cell = {bin, segment, activation}; diff --git a/Tests/UnitTests/Examples/Io/Csv/MeasurementReaderWriterTests.cpp b/Tests/UnitTests/Examples/Io/Csv/MeasurementReaderWriterTests.cpp index 7c09c027ccb..a895e61b64c 100644 --- a/Tests/UnitTests/Examples/Io/Csv/MeasurementReaderWriterTests.cpp +++ b/Tests/UnitTests/Examples/Io/Csv/MeasurementReaderWriterTests.cpp @@ -60,8 +60,8 @@ BOOST_AUTO_TEST_CASE(CsvMeasurementRoundTrip) { ActsExamples::Cluster cl; - using Bin2D = ActsFatras::Channelizer::Bin2D; - using Seg2D = ActsFatras::Channelizer::Segment2D; + using Bin2D = ActsFatras::Segmentizer::Bin2D; + using Seg2D = ActsFatras::Segmentizer::Segment2D; // We have two cluster shapes which are displaced randomly const auto o = disti(gen); diff --git a/Tests/UnitTests/Fatras/Digitization/CMakeLists.txt b/Tests/UnitTests/Fatras/Digitization/CMakeLists.txt index fe1855cb233..8737e163dfc 100644 --- a/Tests/UnitTests/Fatras/Digitization/CMakeLists.txt +++ b/Tests/UnitTests/Fatras/Digitization/CMakeLists.txt @@ -5,3 +5,4 @@ add_unittest(FatrasChannelizer ChannelizerTests.cpp) add_unittest(FatrasPlanarSurfaceDrift PlanarSurfaceDriftTests.cpp) add_unittest(FatrasPlanarSurfaceMask PlanarSurfaceMaskTests.cpp) add_unittest(FatrasUncorrelatedHitSmearer UncorrelatedHitSmearerTests.cpp) +add_unittest(FatrasSegmentizer SegmentizerTests.cpp) diff --git a/Tests/UnitTests/Fatras/Digitization/ChannelizerTests.cpp b/Tests/UnitTests/Fatras/Digitization/ChannelizerTests.cpp index 5c4408f43e1..28748dc1804 100644 --- a/Tests/UnitTests/Fatras/Digitization/ChannelizerTests.cpp +++ b/Tests/UnitTests/Fatras/Digitization/ChannelizerTests.cpp @@ -1,218 +1,143 @@ // This file is part of the Acts project. // -// Copyright (C) 2020 CERN for the benefit of the Acts project +// Copyright (C) 2023 CERN for the benefit of the Acts project // // This Source Code Form is subject to the terms of the Mozilla Public // License, v. 2.0. If a copy of the MPL was not distributed with this // file, You can obtain one at http://mozilla.org/MPL/2.0/. -#include #include -#include "Acts/Definitions/Algebra.hpp" -#include "Acts/Geometry/GeometryContext.hpp" -#include "Acts/Surfaces/DiscBounds.hpp" -#include "Acts/Surfaces/DiscSurface.hpp" -#include "Acts/Surfaces/PlanarBounds.hpp" +#include "Acts/Definitions/Units.hpp" #include "Acts/Surfaces/PlaneSurface.hpp" -#include "Acts/Surfaces/RadialBounds.hpp" -#include "Acts/Surfaces/RectangleBounds.hpp" -#include "Acts/Surfaces/Surface.hpp" #include "Acts/Utilities/BinUtility.hpp" -#include "Acts/Utilities/BinningType.hpp" #include "ActsFatras/Digitization/Channelizer.hpp" +#include "ActsFatras/Digitization/PlanarSurfaceDrift.hpp" +#include "ActsFatras/Digitization/PlanarSurfaceMask.hpp" -#include -#include -#include -#include -#include -#include -#include - -#include "DigitizationCsvOutput.hpp" -#include "PlanarSurfaceTestBeds.hpp" - -namespace bdata = boost::unit_test::data; - -namespace ActsFatras { - -BOOST_AUTO_TEST_SUITE(Digitization) - -BOOST_AUTO_TEST_CASE(ChannelizerCartesian) { - Acts::GeometryContext geoCtx; - - auto rectangleBounds = std::make_shared(1., 1.); - auto planeSurface = Acts::Surface::makeShared( - Acts::Transform3::Identity(), rectangleBounds); - - // The segmentation - Acts::BinUtility pixelated(20, -1., 1., Acts::open, Acts::binX); - pixelated += Acts::BinUtility(20, -1., 1., Acts::open, Acts::binY); - - Channelizer cl; - - // Test: Normal hit into the surface - Acts::Vector2 nPosition(0.37, 0.76); - auto nSegments = - cl.segments(geoCtx, *planeSurface, pixelated, {nPosition, nPosition}); - BOOST_CHECK(nSegments.size() == 1); - BOOST_CHECK(nSegments[0].bin[0] == 13); - BOOST_CHECK(nSegments[0].bin[1] == 17); - - // Test: Inclined hit into the surface - negative x direction - Acts::Vector2 ixPositionS(0.37, 0.76); - Acts::Vector2 ixPositionE(0.02, 0.73); - auto ixSegments = - cl.segments(geoCtx, *planeSurface, pixelated, {ixPositionS, ixPositionE}); - BOOST_CHECK(ixSegments.size() == 4); - - // Test: Inclined hit into the surface - positive y direction - Acts::Vector2 iyPositionS(0.37, 0.76); - Acts::Vector2 iyPositionE(0.39, 0.91); - auto iySegments = - cl.segments(geoCtx, *planeSurface, pixelated, {iyPositionS, iyPositionE}); - BOOST_CHECK(iySegments.size() == 3); - - // Test: Inclined hit into the surface - x/y direction - Acts::Vector2 ixyPositionS(-0.27, 0.76); - Acts::Vector2 ixyPositionE(-0.02, -0.73); - auto ixySegments = cl.segments(geoCtx, *planeSurface, pixelated, - {ixyPositionS, ixyPositionE}); - BOOST_CHECK(ixySegments.size() == 18); +#include + +using namespace Acts::UnitLiterals; + +struct Helper { + std::shared_ptr surface; + Acts::BinUtility segmentation; + + Acts::GeometryContext gctx{}; + double thickness = 125_um; + Acts::Vector3 driftDir = Acts::Vector3::Zero(); + + ActsFatras::Channelizer channelizer; + + Helper() { + surface = Acts::Surface::makeShared( + Acts::Vector3::Zero(), Acts::Vector3{0.0, 0.0, 1.0}); + + float pitchSize = 50_um; + float min = -200_um; + float max = 200_um; + int bins = (max - min) / pitchSize; + segmentation = Acts::BinUtility(bins, min, max, Acts::BinningOption::open, + Acts::BinningValue::binX); + segmentation += Acts::BinUtility(bins, min, max, Acts::BinningOption::open, + Acts::BinningValue::binY); + } + + auto channelize(const Acts::Vector3 &pos3, const Acts::Vector3 &dir3) const { + Acts::Vector4 pos4 = Acts::Vector4::Zero(); + pos4.segment<3>(1) = pos3; + Acts::Vector4 mom4 = Acts::Vector4::Zero(); + mom4.segment<3>(1) = dir3; + ActsFatras::Hit hit({}, {}, pos4, mom4, mom4); + auto res = channelizer.channelize(hit, *surface, gctx, driftDir, + segmentation, thickness); + BOOST_REQUIRE(res.ok()); + return *res; + } +}; + +BOOST_AUTO_TEST_CASE(test_upright_particle) { + Helper helper; + + Acts::Vector3 pos3 = Acts::Vector3{10_um, 10_um, 0.0}; + Acts::Vector3 dir3 = Acts::Vector3{0.0, 0.0, helper.thickness}.normalized(); + + auto segments = helper.channelize(pos3, dir3); + + BOOST_CHECK(segments.size() == 1); + BOOST_CHECK_CLOSE(segments[0].activation, helper.thickness, 1.e-8); } -BOOST_AUTO_TEST_CASE(ChannelizerPolarRadial) { - Acts::GeometryContext geoCtx; - - auto radialBounds = - std::make_shared(5., 10., 0.25, 0.); - auto radialDisc = Acts::Surface::makeShared( - Acts::Transform3::Identity(), radialBounds); - - // The segmentation - Acts::BinUtility strips(2, 5., 10., Acts::open, Acts::binR); - strips += Acts::BinUtility(250, -0.25, 0.25, Acts::open, Acts::binPhi); - - Channelizer cl; - - // Test: Normal hit into the surface - Acts::Vector2 nPosition(6.76, 0.5); - auto nSegments = - cl.segments(geoCtx, *radialDisc, strips, {nPosition, nPosition}); - BOOST_CHECK(nSegments.size() == 1); - BOOST_CHECK(nSegments[0].bin[0] == 0); - BOOST_CHECK(nSegments[0].bin[1] == 161); - - // Test: now opver more phi strips - Acts::Vector2 sPositionS(6.76, 0.5); - Acts::Vector2 sPositionE(7.03, -0.3); - auto sSegment = - cl.segments(geoCtx, *radialDisc, strips, {sPositionS, sPositionE}); - BOOST_CHECK(sSegment.size() == 59); - - // Test: jump over R boundary, but stay in phi bin - sPositionS = Acts::Vector2(6.76, 0.); - sPositionE = Acts::Vector2(7.83, 0.); - sSegment = cl.segments(geoCtx, *radialDisc, strips, {sPositionS, sPositionE}); - BOOST_CHECK(sSegment.size() == 2); +BOOST_AUTO_TEST_CASE(test_tilted_particle) { + Helper helper; + + const double disp = 10_um; + + Acts::Vector3 hitPosition = Acts::Vector3{10_um, 10_um, 0.0}; + Acts::Vector3 hitDirection = + Acts::Vector3({disp, 0.0, helper.thickness}).normalized(); + + auto segments = helper.channelize(hitPosition, hitDirection); + + BOOST_CHECK(segments.size() == 1); + BOOST_CHECK_CLOSE(segments[0].activation, std::hypot(disp, helper.thickness), + 1.e-8); } -/// Unit test for testing the Channelizer -BOOST_DATA_TEST_CASE(RandomChannelizerTest, - bdata::random(0., 1.) ^ bdata::random(0., 1.) ^ - bdata::random(0., 1.) ^ bdata::random(0., 1.) ^ - bdata::xrange(25), - startR0, startR1, endR0, endR1, index) { - Acts::GeometryContext geoCtx; - Channelizer cl; - - // Test beds with random numbers generated inside - PlanarSurfaceTestBeds pstd; - auto testBeds = pstd(1.); - - DigitizationCsvOutput csvHelper; - - for (const auto& tb : testBeds) { - const auto& name = std::get<0>(tb); - const auto* surface = (std::get<1>(tb)).get(); - const auto& segmentation = std::get<2>(tb); - const auto& randomizer = std::get<3>(tb); - - if (index == 0) { - std::ofstream shape; - std::ofstream grid; - const auto centerXY = surface->center(geoCtx).segment<2>(0); - // 0 - write the shape - shape.open("Channelizer" + name + "Borders.csv"); - if (surface->type() == Acts::Surface::Plane) { - const auto* pBounds = - static_cast(&(surface->bounds())); - csvHelper.writePolygon(shape, pBounds->vertices(1), -centerXY); - } else if (surface->type() == Acts::Surface::Disc) { - const auto* dBounds = - static_cast(&(surface->bounds())); - csvHelper.writePolygon(shape, dBounds->vertices(72), -centerXY); - } - // 1 - write the grid - grid.open("Channelizer" + name + "Grid.csv"); - if (segmentation.binningData()[0].binvalue == Acts::binX && - segmentation.binningData()[1].binvalue == Acts::binY) { - double bxmin = segmentation.binningData()[0].min; - double bxmax = segmentation.binningData()[0].max; - double bymin = segmentation.binningData()[1].min; - double bymax = segmentation.binningData()[1].max; - const auto& xboundaries = segmentation.binningData()[0].boundaries(); - const auto& yboundaries = segmentation.binningData()[1].boundaries(); - for (const auto xval : xboundaries) { - csvHelper.writeLine(grid, {xval, bymin}, {xval, bymax}); - } - for (const auto yval : yboundaries) { - csvHelper.writeLine(grid, {bxmin, yval}, {bxmax, yval}); - } - } else if (segmentation.binningData()[0].binvalue == Acts::binR && - segmentation.binningData()[1].binvalue == Acts::binPhi) { - double brmin = segmentation.binningData()[0].min; - double brmax = segmentation.binningData()[0].max; - double bphimin = segmentation.binningData()[1].min; - double bphimax = segmentation.binningData()[1].max; - const auto& rboundaries = segmentation.binningData()[0].boundaries(); - const auto& phiboundaries = segmentation.binningData()[1].boundaries(); - for (const auto r : rboundaries) { - csvHelper.writeArc(grid, r, bphimin, bphimax); - } - for (const auto phi : phiboundaries) { - double cphi = std::cos(phi); - double sphi = std::sin(phi); - csvHelper.writeLine(grid, {brmin * cphi, brmin * sphi}, - {brmax * cphi, brmax * sphi}); - } - } - } - - auto start = randomizer(startR0, startR1); - auto end = randomizer(endR0, endR1); - - std::ofstream segments; - segments.open("Channelizer" + name + "Segments_n" + std::to_string(index) + - ".csv"); - - std::ofstream cluster; - cluster.open("Channelizer" + name + "Cluster_n" + std::to_string(index) + - ".csv"); - - /// Run the channelizer - auto cSegement = cl.segments(geoCtx, *surface, segmentation, {start, end}); - - for (const auto& cs : cSegement) { - csvHelper.writeLine(segments, cs.path2D[0], cs.path2D[1]); - } - - segments.close(); - cluster.close(); - } +BOOST_AUTO_TEST_CASE(test_more_tilted_particle) { + Helper helper; + + const double disp = 50_um; + + Acts::Vector3 hitPosition = Acts::Vector3{10_um, 10_um, 0.0}; + Acts::Vector3 hitDirection = + Acts::Vector3{disp, 0.0, helper.thickness}.normalized(); + + auto segments = helper.channelize(hitPosition, hitDirection); + + BOOST_CHECK(segments.size() == 2); + auto sum = + std::accumulate(segments.begin(), segments.end(), 0.0, + [](double s, auto seg) { return s + seg.activation; }); + BOOST_CHECK_CLOSE(sum, std::hypot(disp, helper.thickness), 1.e-8); +} + +// This should go directly up on the segment border +BOOST_AUTO_TEST_CASE(test_pathological_upright_particle) { + Helper helper; + + Acts::Vector3 hitPosition = Acts::Vector3{0.0, 10_um, 0.0}; + Acts::Vector3 hitDirection = + Acts::Vector3{0.0, 0.0, helper.thickness}.normalized(); + + auto segments = helper.channelize(hitPosition, hitDirection); + + BOOST_CHECK(segments.size() == 1); + BOOST_CHECK_CLOSE(segments[0].activation, helper.thickness, 1.e-8); } -BOOST_AUTO_TEST_SUITE_END() +// This should go directly up on the segment border +// TODO why does this does not activate both cells with half of the path??? +BOOST_AUTO_TEST_CASE(test_pathological_tilted_particle) { + Helper helper; + + double disp = 2.0_um; + + Acts::Vector3 hitPosition = Acts::Vector3{-0.5 * disp, 10_um, 0.0}; + Acts::Vector3 hitDirection = + Acts::Vector3{disp, 0.0, helper.thickness}.normalized(); -} // namespace ActsFatras + auto segments = helper.channelize(hitPosition, hitDirection); + + std::cout << "Segments:\n"; + for (const auto &seg : segments) { + std::cout << " - (" << seg.bin[0] << ", " << seg.bin[1] + << "), activation: " << seg.activation << "\n"; + } + + BOOST_CHECK(segments.size() == 2); + auto sum = + std::accumulate(segments.begin(), segments.end(), 0.0, + [](double s, auto seg) { return s + seg.activation; }); + BOOST_CHECK_CLOSE(sum, std::hypot(disp, helper.thickness), 1.e-8); +} diff --git a/Tests/UnitTests/Fatras/Digitization/SegmentizerTests.cpp b/Tests/UnitTests/Fatras/Digitization/SegmentizerTests.cpp new file mode 100644 index 00000000000..e8013f538b3 --- /dev/null +++ b/Tests/UnitTests/Fatras/Digitization/SegmentizerTests.cpp @@ -0,0 +1,218 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2020 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include +#include + +#include "Acts/Definitions/Algebra.hpp" +#include "Acts/Geometry/GeometryContext.hpp" +#include "Acts/Surfaces/DiscBounds.hpp" +#include "Acts/Surfaces/DiscSurface.hpp" +#include "Acts/Surfaces/PlanarBounds.hpp" +#include "Acts/Surfaces/PlaneSurface.hpp" +#include "Acts/Surfaces/RadialBounds.hpp" +#include "Acts/Surfaces/RectangleBounds.hpp" +#include "Acts/Surfaces/Surface.hpp" +#include "Acts/Utilities/BinUtility.hpp" +#include "Acts/Utilities/BinningType.hpp" +#include "ActsFatras/Digitization/Segmentizer.hpp" + +#include +#include +#include +#include +#include +#include +#include + +#include "DigitizationCsvOutput.hpp" +#include "PlanarSurfaceTestBeds.hpp" + +namespace bdata = boost::unit_test::data; + +namespace ActsFatras { + +BOOST_AUTO_TEST_SUITE(Digitization) + +BOOST_AUTO_TEST_CASE(SegmentizerCartesian) { + Acts::GeometryContext geoCtx; + + auto rectangleBounds = std::make_shared(1., 1.); + auto planeSurface = Acts::Surface::makeShared( + Acts::Transform3::Identity(), rectangleBounds); + + // The segmentation + Acts::BinUtility pixelated(20, -1., 1., Acts::open, Acts::binX); + pixelated += Acts::BinUtility(20, -1., 1., Acts::open, Acts::binY); + + Segmentizer cl; + + // Test: Normal hit into the surface + Acts::Vector2 nPosition(0.37, 0.76); + auto nSegments = + cl.segments(geoCtx, *planeSurface, pixelated, {nPosition, nPosition}); + BOOST_CHECK(nSegments.size() == 1); + BOOST_CHECK(nSegments[0].bin[0] == 13); + BOOST_CHECK(nSegments[0].bin[1] == 17); + + // Test: Inclined hit into the surface - negative x direction + Acts::Vector2 ixPositionS(0.37, 0.76); + Acts::Vector2 ixPositionE(0.02, 0.73); + auto ixSegments = + cl.segments(geoCtx, *planeSurface, pixelated, {ixPositionS, ixPositionE}); + BOOST_CHECK(ixSegments.size() == 4); + + // Test: Inclined hit into the surface - positive y direction + Acts::Vector2 iyPositionS(0.37, 0.76); + Acts::Vector2 iyPositionE(0.39, 0.91); + auto iySegments = + cl.segments(geoCtx, *planeSurface, pixelated, {iyPositionS, iyPositionE}); + BOOST_CHECK(iySegments.size() == 3); + + // Test: Inclined hit into the surface - x/y direction + Acts::Vector2 ixyPositionS(-0.27, 0.76); + Acts::Vector2 ixyPositionE(-0.02, -0.73); + auto ixySegments = cl.segments(geoCtx, *planeSurface, pixelated, + {ixyPositionS, ixyPositionE}); + BOOST_CHECK(ixySegments.size() == 18); +} + +BOOST_AUTO_TEST_CASE(SegmentizerPolarRadial) { + Acts::GeometryContext geoCtx; + + auto radialBounds = + std::make_shared(5., 10., 0.25, 0.); + auto radialDisc = Acts::Surface::makeShared( + Acts::Transform3::Identity(), radialBounds); + + // The segmentation + Acts::BinUtility strips(2, 5., 10., Acts::open, Acts::binR); + strips += Acts::BinUtility(250, -0.25, 0.25, Acts::open, Acts::binPhi); + + Segmentizer cl; + + // Test: Normal hit into the surface + Acts::Vector2 nPosition(6.76, 0.5); + auto nSegments = + cl.segments(geoCtx, *radialDisc, strips, {nPosition, nPosition}); + BOOST_CHECK(nSegments.size() == 1); + BOOST_CHECK(nSegments[0].bin[0] == 0); + BOOST_CHECK(nSegments[0].bin[1] == 161); + + // Test: now opver more phi strips + Acts::Vector2 sPositionS(6.76, 0.5); + Acts::Vector2 sPositionE(7.03, -0.3); + auto sSegment = + cl.segments(geoCtx, *radialDisc, strips, {sPositionS, sPositionE}); + BOOST_CHECK(sSegment.size() == 59); + + // Test: jump over R boundary, but stay in phi bin + sPositionS = Acts::Vector2(6.76, 0.); + sPositionE = Acts::Vector2(7.83, 0.); + sSegment = cl.segments(geoCtx, *radialDisc, strips, {sPositionS, sPositionE}); + BOOST_CHECK(sSegment.size() == 2); +} + +/// Unit test for testing the Segmentizer +BOOST_DATA_TEST_CASE(RandomSegmentizerTest, + bdata::random(0., 1.) ^ bdata::random(0., 1.) ^ + bdata::random(0., 1.) ^ bdata::random(0., 1.) ^ + bdata::xrange(25), + startR0, startR1, endR0, endR1, index) { + Acts::GeometryContext geoCtx; + Segmentizer cl; + + // Test beds with random numbers generated inside + PlanarSurfaceTestBeds pstd; + auto testBeds = pstd(1.); + + DigitizationCsvOutput csvHelper; + + for (const auto& tb : testBeds) { + const auto& name = std::get<0>(tb); + const auto* surface = (std::get<1>(tb)).get(); + const auto& segmentation = std::get<2>(tb); + const auto& randomizer = std::get<3>(tb); + + if (index == 0) { + std::ofstream shape; + std::ofstream grid; + const auto centerXY = surface->center(geoCtx).segment<2>(0); + // 0 - write the shape + shape.open("Segmentizer" + name + "Borders.csv"); + if (surface->type() == Acts::Surface::Plane) { + const auto* pBounds = + static_cast(&(surface->bounds())); + csvHelper.writePolygon(shape, pBounds->vertices(1), -centerXY); + } else if (surface->type() == Acts::Surface::Disc) { + const auto* dBounds = + static_cast(&(surface->bounds())); + csvHelper.writePolygon(shape, dBounds->vertices(72), -centerXY); + } + // 1 - write the grid + grid.open("Segmentizer" + name + "Grid.csv"); + if (segmentation.binningData()[0].binvalue == Acts::binX && + segmentation.binningData()[1].binvalue == Acts::binY) { + double bxmin = segmentation.binningData()[0].min; + double bxmax = segmentation.binningData()[0].max; + double bymin = segmentation.binningData()[1].min; + double bymax = segmentation.binningData()[1].max; + const auto& xboundaries = segmentation.binningData()[0].boundaries(); + const auto& yboundaries = segmentation.binningData()[1].boundaries(); + for (const auto xval : xboundaries) { + csvHelper.writeLine(grid, {xval, bymin}, {xval, bymax}); + } + for (const auto yval : yboundaries) { + csvHelper.writeLine(grid, {bxmin, yval}, {bxmax, yval}); + } + } else if (segmentation.binningData()[0].binvalue == Acts::binR && + segmentation.binningData()[1].binvalue == Acts::binPhi) { + double brmin = segmentation.binningData()[0].min; + double brmax = segmentation.binningData()[0].max; + double bphimin = segmentation.binningData()[1].min; + double bphimax = segmentation.binningData()[1].max; + const auto& rboundaries = segmentation.binningData()[0].boundaries(); + const auto& phiboundaries = segmentation.binningData()[1].boundaries(); + for (const auto r : rboundaries) { + csvHelper.writeArc(grid, r, bphimin, bphimax); + } + for (const auto phi : phiboundaries) { + double cphi = std::cos(phi); + double sphi = std::sin(phi); + csvHelper.writeLine(grid, {brmin * cphi, brmin * sphi}, + {brmax * cphi, brmax * sphi}); + } + } + } + + auto start = randomizer(startR0, startR1); + auto end = randomizer(endR0, endR1); + + std::ofstream segments; + segments.open("Segmentizer" + name + "Segments_n" + std::to_string(index) + + ".csv"); + + std::ofstream cluster; + cluster.open("Segmentizer" + name + "Cluster_n" + std::to_string(index) + + ".csv"); + + /// Run the Segmentizer + auto cSegement = cl.segments(geoCtx, *surface, segmentation, {start, end}); + + for (const auto& cs : cSegement) { + csvHelper.writeLine(segments, cs.path2D[0], cs.path2D[1]); + } + + segments.close(); + cluster.close(); + } +} + +BOOST_AUTO_TEST_SUITE_END() + +} // namespace ActsFatras diff --git a/Tests/UnitTests/Plugins/ExaTrkX/ExaTrkXBoostTrackBuildingTests.cpp b/Tests/UnitTests/Plugins/ExaTrkX/ExaTrkXBoostTrackBuildingTests.cpp index 1adddf15638..b5ca404006a 100644 --- a/Tests/UnitTests/Plugins/ExaTrkX/ExaTrkXBoostTrackBuildingTests.cpp +++ b/Tests/UnitTests/Plugins/ExaTrkX/ExaTrkXBoostTrackBuildingTests.cpp @@ -9,16 +9,21 @@ #include #include "Acts/Plugins/ExaTrkX/BoostTrackBuilding.hpp" +#include "Acts/Plugins/ExaTrkX/detail/BoostTrackBuildingUtils.hpp" #include "Acts/Plugins/ExaTrkX/detail/TensorVectorConversion.hpp" #include +#include + BOOST_AUTO_TEST_CASE(test_track_building) { // Make some spacepoint IDs // The spacepoint ids are [100, 101, 102, ...] // They should not be zero based to check if the thing also works if the // spacepoint IDs do not match the node IDs used for the edges std::vector spacepointIds(16); + auto nodes = torch::rand({16, 3}); + std::iota(spacepointIds.begin(), spacepointIds.end(), 100); // Build 4 tracks with 4 hits @@ -46,7 +51,8 @@ BOOST_AUTO_TEST_CASE(test_track_building) { auto logger = Acts::getDefaultLogger("TestLogger", Acts::Logging::ERROR); Acts::BoostTrackBuilding trackBuilder(std::move(logger)); - auto testTracks = trackBuilder({}, edgeTensor, dummyWeights, spacepointIds); + auto testTracks = + trackBuilder(nodes, edgeTensor, dummyWeights, spacepointIds); // Sort tracks, so we can find them std::for_each(testTracks.begin(), testTracks.end(), @@ -60,3 +66,79 @@ BOOST_AUTO_TEST_CASE(test_track_building) { BOOST_CHECK(found != testTracks.end()); } } + +struct EdgeProperty { + float weight; +}; + +BOOST_AUTO_TEST_CASE(test_graph_cleaning_no_cleaning) { + using Graph = + boost::adjacency_list; + + Graph graph; + + // Add one clean track + boost::add_edge(0, 1, EdgeProperty{1.0}, graph); + boost::add_edge(1, 2, EdgeProperty{1.0}, graph); + boost::add_edge(2, 3, EdgeProperty{1.0}, graph); + boost::add_edge(3, 4, EdgeProperty{1.0}, graph); + + // Add another clean track, but with one weak edge (should have no effect) + boost::add_edge(5, 6, EdgeProperty{1.0}, graph); + boost::add_edge(6, 7, EdgeProperty{0.5}, graph); + boost::add_edge(7, 8, EdgeProperty{1.0}, graph); + boost::add_edge(8, 9, EdgeProperty{1.0}, graph); + + std::size_t numEdgesBefore = boost::num_edges(graph); + std::vector c(boost::num_vertices(graph)); + std::size_t numConnectedBefore = boost::connected_components(graph, c.data()); + + Acts::detail::cleanSubgraphs(graph); + + BOOST_CHECK(numEdgesBefore == boost::num_edges(graph)); + BOOST_CHECK(numConnectedBefore == + boost::connected_components(graph, c.data())); +} + +BOOST_AUTO_TEST_CASE(test_graph_cleaning_one_branch) { + using Graph = + boost::adjacency_list; + + Graph graph; + + // Add one clean track + boost::add_edge(0, 1, EdgeProperty{1.0}, graph); + boost::add_edge(1, 2, EdgeProperty{1.0}, graph); + boost::add_edge(2, 3, EdgeProperty{1.0}, graph); + boost::add_edge(3, 4, EdgeProperty{1.0}, graph); + + // Add another clean track + // Should be branched in 5-9 and 10-12 + boost::add_edge(5, 6, EdgeProperty{1.0}, graph); + + boost::add_edge(6, 7, EdgeProperty{1.0}, graph); + boost::add_edge(7, 8, EdgeProperty{1.0}, graph); + boost::add_edge(8, 9, EdgeProperty{1.0}, graph); + + boost::add_edge(6, 10, EdgeProperty{0.5}, graph); + boost::add_edge(10, 11, EdgeProperty{1.0}, graph); + boost::add_edge(11, 12, EdgeProperty{1.0}, graph); + + std::size_t numEdgesBefore = boost::num_edges(graph); + std::vector c(boost::num_vertices(graph)); + std::size_t numConnectedBefore = boost::connected_components(graph, c.data()); + + auto logger = Acts::getDefaultLogger("TestLogger", Acts::Logging::VERBOSE); + Acts::detail::cleanSubgraphs(graph, *logger); + + std::cout << "edges count " << numEdgesBefore << " -> " + << boost::num_edges(graph) << std::endl; + std::cout << "connected cmp count " << numConnectedBefore << " -> " + << boost::connected_components(graph, c.data()) << std::endl; + + BOOST_CHECK(boost::num_edges(graph) == numEdgesBefore - 1); + BOOST_CHECK(boost::connected_components(graph, c.data()) == + numConnectedBefore + 1); +} diff --git a/Tests/UnitTests/Plugins/ExaTrkX/ExaTrkXEdgeBuildingTests.cpp b/Tests/UnitTests/Plugins/ExaTrkX/ExaTrkXEdgeBuildingTests.cpp new file mode 100644 index 00000000000..01068c7c4ad --- /dev/null +++ b/Tests/UnitTests/Plugins/ExaTrkX/ExaTrkXEdgeBuildingTests.cpp @@ -0,0 +1,257 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2022 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include + +#include "Acts/Plugins/ExaTrkX/detail/CantorEdge.hpp" +#include "Acts/Plugins/ExaTrkX/detail/TensorVectorConversion.hpp" +#include "Acts/Plugins/ExaTrkX/detail/buildEdges.hpp" +#include "Acts/Utilities/ContainerPrinter.hpp" + +#include +#include + +#include +#include + +using CantorPair = Acts::detail::CantorEdge; + +#define PRINT 0 + +float distance(const at::Tensor &a, const at::Tensor &b) { + assert(a.sizes() == b.sizes()); + assert(a.sizes().size() == 1); + + return std::sqrt(((a - b) * (a - b)).sum().item().to()); +} + +#if PRINT +std::ostream &operator<<(std::ostream &os, CantorPair p) { + auto [a, b] = p.inverse(); + os << "(" << a << "," << b << ")"; + return os; +} +#endif + +template +void test_random_graph(int emb_dim, int n_nodes, float r, int knn, + const edge_builder_t &edgeBuilder) { + // Create a random point cloud + auto random_features = at::randn({n_nodes, emb_dim}); + + // Generate the truth via brute-force + Eigen::MatrixXf distance_matrix(n_nodes, n_nodes); + + std::vector edges_ref_cantor; + std::vector edge_counts(n_nodes, 0); + + for (int i = 0; i < n_nodes; ++i) { + for (int j = i; j < n_nodes; ++j) { + const auto d = distance(random_features[i], random_features[j]); + distance_matrix(i, j) = d; + distance_matrix(j, i) = d; + + if (d < r && i != j) { + edges_ref_cantor.emplace_back(i, j); + edge_counts[i]++; + } + } + } + + const auto max_edges = + *std::max_element(edge_counts.begin(), edge_counts.end()); + + // If this is not the case, the test is ill-formed + // knn specifies how many edges can be found by the function at max. Thus we + // should design the test in a way, that our brute-force test algorithm does + // not find more edges then the algorithm that we test against it can find + BOOST_REQUIRE(max_edges <= knn); + + // Run the edge building + auto edges_test = edgeBuilder(random_features, r, knn); + + // Map the edges to cantor pairs + std::vector edges_test_cantor; + + for (int i = 0; i < edges_test.size(1); ++i) { + const auto a = edges_test[0][i].template item(); + const auto b = edges_test[1][i].template item(); + edges_test_cantor.push_back(a < b ? CantorPair(a, b) : CantorPair(b, a)); + } + + std::sort(edges_ref_cantor.begin(), edges_ref_cantor.end()); + std::sort(edges_test_cantor.begin(), edges_test_cantor.end()); + +#if PRINT + std::cout << "test size " << edges_test_cantor.size() << std::endl; + std::cout << "ref size " << edges_ref_cantor.size() << std::endl; + std::cout << "test: " << Acts::ContainerPrinter(edges_test_cantor, 10) + << std::endl; + std::cout << "ref: " << Acts::ContainerPrinter(edges_ref_cantor, 10) + << std::endl; +#endif + + // Check + BOOST_CHECK(edges_ref_cantor.size() == edges_test_cantor.size()); + BOOST_CHECK(std::equal(edges_test_cantor.begin(), edges_test_cantor.end(), + edges_ref_cantor.begin())); +} + +BOOST_AUTO_TEST_CASE(test_cantor_pair_functions) { + int a = 345; + int b = 23; + const auto [aa, bb] = CantorPair(a, b).inverse(); + BOOST_CHECK(a == aa); + BOOST_CHECK(b == bb); +} + +const int emb_dim = 3; +const int n_nodes = 20; +const float r = 1.5; +const int knn = 50; +const int seed = 42; + +BOOST_AUTO_TEST_CASE(test_random_graph_edge_building_cuda, + *boost::unit_test::precondition([](auto) { + return torch::cuda::is_available(); + })) { + torch::manual_seed(seed); + + auto cudaEdgeBuilder = [](auto &features, auto radius, auto k) { + auto features_cuda = features.to(torch::kCUDA); + return Acts::detail::buildEdgesFRNN(features_cuda, radius, k); + }; + + test_random_graph(emb_dim, n_nodes, r, knn, cudaEdgeBuilder); +} + +BOOST_AUTO_TEST_CASE(test_random_graph_edge_building_kdtree) { + torch::manual_seed(seed); + + auto cpuEdgeBuilder = [](auto &features, auto radius, auto k) { + auto features_cpu = features.to(torch::kCPU); + return Acts::detail::buildEdgesKDTree(features_cpu, radius, k); + }; + + test_random_graph(emb_dim, n_nodes, r, knn, cpuEdgeBuilder); +} + +BOOST_AUTO_TEST_CASE(test_self_loop_removal) { + // clang-format off + std::vector edges = { + 1,1, + 2,3, + 2,2, + 5,4, + }; + // clang-format on + + auto opts = torch::TensorOptions().dtype(torch::kInt64); + const auto edgeTensor = + torch::from_blob(edges.data(), {static_cast(edges.size() / 2), 2}, + opts) + .transpose(0, 1); + + const auto withoutSelfLoops = + Acts::detail::postprocessEdgeTensor(edgeTensor, true, false, false) + .transpose(1, 0) + .flatten(); + + const std::vector postEdges( + withoutSelfLoops.data_ptr(), + withoutSelfLoops.data_ptr() + withoutSelfLoops.numel()); + + // clang-format off + const std::vector ref = { + 2,3, + 5,4, + }; + // clang-format on + + BOOST_CHECK(ref == postEdges); +} + +BOOST_AUTO_TEST_CASE(test_duplicate_removal) { + // clang-format off + std::vector edges = { + 1,2, + 2,1, // duplicate, flipped + 3,2, + 3,2, // duplicate, not flipped + 7,6, // should be flipped + }; + // clang-format on + + auto opts = torch::TensorOptions().dtype(torch::kInt64); + const auto edgeTensor = + torch::from_blob(edges.data(), {static_cast(edges.size() / 2), 2}, + opts) + .transpose(0, 1); + + const auto withoutDups = + Acts::detail::postprocessEdgeTensor(edgeTensor, false, true, false) + .transpose(1, 0) + .flatten(); + + const std::vector postEdges( + withoutDups.data_ptr(), + withoutDups.data_ptr() + withoutDups.numel()); + + // clang-format off + const std::vector ref = { + 1,2, + 2,3, + 6,7, + }; + // clang-format on + + BOOST_CHECK(ref == postEdges); +} + +BOOST_AUTO_TEST_CASE(test_random_flip) { + torch::manual_seed(seed); + + // clang-format off + std::vector edges = { + 1,2, + 2,3, + 3,4, + 4,5, + }; + // clang-format on + + auto opts = torch::TensorOptions().dtype(torch::kInt64); + const auto edgeTensor = + torch::from_blob(edges.data(), {static_cast(edges.size() / 2), 2}, + opts) + .transpose(0, 1); + + const auto flipped = + Acts::detail::postprocessEdgeTensor(edgeTensor, false, false, true) + .transpose(0, 1) + .flatten(); + + const std::vector postEdges( + flipped.data_ptr(), + flipped.data_ptr() + flipped.numel()); + + BOOST_CHECK(postEdges.size() == edges.size()); + for (auto preIt = edges.begin(); preIt != edges.end(); preIt += 2) { + int found = 0; + + for (auto postIt = postEdges.begin(); postIt != postEdges.end(); + postIt += 2) { + bool noflp = (*preIt == *postIt) and *(preIt + 1) == *(postIt + 1); + bool flp = *preIt == *(postIt + 1) and *(preIt + 1) == *(postIt); + + found += (flp or noflp); + } + + BOOST_CHECK(found == 1); + } +} diff --git a/Tests/UnitTests/Plugins/ExaTrkX/ExaTrkXEdgeClassifierTests.cpp b/Tests/UnitTests/Plugins/ExaTrkX/ExaTrkXEdgeClassifierTests.cpp new file mode 100644 index 00000000000..07fccdaf33e --- /dev/null +++ b/Tests/UnitTests/Plugins/ExaTrkX/ExaTrkXEdgeClassifierTests.cpp @@ -0,0 +1,46 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2022 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include + +#include "Acts/Plugins/ExaTrkX/TorchEdgeClassifier.hpp" +#include "Acts/Plugins/ExaTrkX/detail/TensorVectorConversion.hpp" +#include "Acts/Plugins/ExaTrkX/detail/buildEdges.hpp" + +#include +#include + +#include +#include + +using namespace Acts; + +BOOST_AUTO_TEST_CASE(test_model) { + TorchEdgeClassifier::Config cfg; + cfg.modelPath = "/home/benjamin/Desktop/test.pt"; + cfg.cut = 0.0; + cfg.nChunks = 1; + cfg.undirected = false; + + auto logger = Acts::getDefaultLogger("test", Logging::INFO); + TorchEdgeClassifier classifier(cfg, std::move(logger)); + + auto nodes = torch::rand({20, 3}).to(torch::kFloat32); + auto edges = torch::randint(0, 20, {2, 10}); + + auto [n, e, w] = classifier(nodes, edges); + std::cout << std::get<0>(torch::sort(std::any_cast(w))) + << std::endl; + + auto nodes2 = torch::roll(nodes, 1, 0); + auto edges2 = (edges + 1) % 20; + + auto [n2, e2, w2] = classifier(nodes2, edges2); + std::cout << std::get<0>(torch::sort(std::any_cast(w2))) + << std::endl; +} diff --git a/thirdparty/FRNN/CMakeLists.txt b/thirdparty/FRNN/CMakeLists.txt index 9aa70a450aa..a905daf8498 100644 --- a/thirdparty/FRNN/CMakeLists.txt +++ b/thirdparty/FRNN/CMakeLists.txt @@ -18,9 +18,15 @@ set(ACTS_FRNN_GIT_TAG "3e370d8d9073d4e130363faf87d2370598b5fbf2" mark_as_advanced(ACTS_FRNN_GIT_REPOSITORY ACTS_FRNN_GIT_TAG) +set(PATCH ${CMAKE_CURRENT_SOURCE_DIR}/thread_local.patch) + +# Here we apply a patch to add the thread_local keyword to some global variables. +# This enables the use of the library in multi-threaded contexts (I hope). +# Either apply the patch or check if the patch is already applied FetchContent_Declare(frnncontent GIT_REPOSITORY "${ACTS_FRNN_GIT_REPOSITORY}" - GIT_TAG "${ACTS_FRNN_GIT_TAG}" ) + GIT_TAG "${ACTS_FRNN_GIT_TAG}" + PATCH_COMMAND git apply ${PATCH} || git apply --reverse --check ${PATCH}) # FRNN does not provide a CMakeLists.txt, so we use a custom one. Because of this, # we have to implement the populate step manually diff --git a/thirdparty/FRNN/thread_local.patch b/thirdparty/FRNN/thread_local.patch new file mode 100644 index 00000000000..4c76669bd1f --- /dev/null +++ b/thirdparty/FRNN/thread_local.patch @@ -0,0 +1,17 @@ +diff --git a/frnn/csrc/grid/prefix_sum.cu b/frnn/csrc/grid/prefix_sum.cu +index d67ced2..7d332ee 100644 +--- a/frnn/csrc/grid/prefix_sum.cu ++++ b/frnn/csrc/grid/prefix_sum.cu +@@ -161,9 +161,9 @@ inline int floorPow2(int n) { + + #define BLOCK_SIZE 256 + +-int **g_scanBlockSumsInt = 0; +-unsigned int g_numEltsAllocated = 0; +-unsigned int g_numLevelsAllocated = 0; ++thread_local int **g_scanBlockSumsInt = 0; ++thread_local unsigned int g_numEltsAllocated = 0; ++thread_local unsigned int g_numLevelsAllocated = 0; + + bool cudaCheck(cudaError_t status, const std::string &msg) { + if (status != cudaSuccess) {