Skip to content

Commit

Permalink
Merge branch 'main' into intersection-material-assigner-empty-fix
Browse files Browse the repository at this point in the history
  • Loading branch information
AJPfleger authored Apr 30, 2024
2 parents 1ec9819 + 41c0e87 commit 3e6442c
Show file tree
Hide file tree
Showing 11 changed files with 315 additions and 25 deletions.
2 changes: 1 addition & 1 deletion Core/include/Acts/Propagator/EigenStepper.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -431,7 +431,7 @@ class EigenStepper {
std::shared_ptr<const MagneticFieldProvider> m_bField;

/// Overstep limit
double m_overstepLimit;
double m_overstepLimit{};
};

template <typename navigator_t>
Expand Down
51 changes: 36 additions & 15 deletions Core/include/Acts/Propagator/EigenStepper.ipp
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,9 @@
#include "Acts/EventData/TransformationHelpers.hpp"
#include "Acts/Propagator/ConstrainedStep.hpp"
#include "Acts/Propagator/detail/CovarianceEngine.hpp"
#include "Acts/Utilities/QuickMath.hpp"

#include <limits>

template <typename E, typename A>
Acts::EigenStepper<E, A>::EigenStepper(
Expand Down Expand Up @@ -153,7 +156,8 @@ Acts::Result<double> Acts::EigenStepper<E, A>::step(
propagator_state_t& state, const navigator_t& navigator) const {
// Runge-Kutta integrator state
auto& sd = state.stepping.stepData;
double error_estimate = 0.;

double errorEstimate = 0.;
double h2 = 0, half_h = 0;

auto pos = position(state.stepping);
Expand All @@ -172,6 +176,31 @@ Acts::Result<double> Acts::EigenStepper<E, A>::step(
return 0.;
}

const auto calcStepSizeScaling = [&](const double errorEstimate_) -> double {
// For details about these values see ATL-SOFT-PUB-2009-001 for details
constexpr double lower = 0.25;
constexpr double upper = 4.0;
// This is given by the order of the Runge-Kutta method
constexpr double exponent = 0.25;

// Whether to use fast power function if available
constexpr bool tryUseFastPow{false};

double x = state.options.stepTolerance / errorEstimate_;

if constexpr (exponent == 0.25 && !tryUseFastPow) {
// This is 3x faster than std::pow
x = std::sqrt(std::sqrt(x));
} else if constexpr (std::numeric_limits<double>::is_iec559 &&
tryUseFastPow) {
x = fastPow(x, exponent);
} else {
x = std::pow(x, exponent);
}

return std::clamp(x, lower, upper);
};

// The following functor starts to perform a Runge-Kutta step of a certain
// size, going up to the point where it can return an estimate of the local
// integration error. The results are stated in the local variables above,
Expand Down Expand Up @@ -217,12 +246,13 @@ Acts::Result<double> Acts::EigenStepper<E, A>::step(
}

// Compute and check the local integration error estimate
error_estimate =
errorEstimate =
h2 * ((sd.k1 - sd.k2 - sd.k3 + sd.k4).template lpNorm<1>() +
std::abs(sd.kQoP[0] - sd.kQoP[1] - sd.kQoP[2] + sd.kQoP[3]));
error_estimate = std::max(error_estimate, 1e-20);
// Protect against division by zero
errorEstimate = std::max(1e-20, errorEstimate);

return success(error_estimate <= state.options.stepTolerance);
return success(errorEstimate <= state.options.stepTolerance);
};

const double initialH =
Expand All @@ -240,12 +270,7 @@ Acts::Result<double> Acts::EigenStepper<E, A>::step(
break;
}

// double std::sqrt is 3x faster than std::pow
const double stepSizeScaling =
std::clamp(std::sqrt(std::sqrt(state.options.stepTolerance /
std::abs(2. * error_estimate))),
0.25, 4.0);
h *= stepSizeScaling;
h *= calcStepSizeScaling(2 * errorEstimate);

// If step size becomes too small the particle remains at the initial
// place
Expand Down Expand Up @@ -321,11 +346,7 @@ Acts::Result<double> Acts::EigenStepper<E, A>::step(
state.stepping.derivative.template segment<3>(4) = sd.k4;
}
state.stepping.pathAccumulated += h;
// double std::sqrt is 3x faster than std::pow
const double stepSizeScaling =
std::clamp(std::sqrt(std::sqrt(state.options.stepTolerance /
std::abs(error_estimate))),
0.25, 4.0);
const double stepSizeScaling = calcStepSizeScaling(errorEstimate);
const double nextAccuracy = std::abs(h * stepSizeScaling);
const double previousAccuracy = std::abs(state.stepping.stepSize.accuracy());
const double initialStepLength = std::abs(initialH);
Expand Down
90 changes: 90 additions & 0 deletions Core/include/Acts/Utilities/QuickMath.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,90 @@
// This file is part of the Acts project.
//
// Copyright (C) 2024 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 <cmath>
#include <cstdint>
#include <limits>
#include <type_traits>

namespace Acts {

/// @brief Fast inverse square root function
/// Taken from https://en.wikipedia.org/wiki/Fast_inverse_square_root
/// @param x the number to calculate the inverse square root of
/// @param iterations the number of newton iterations to perform
template <typename T>
constexpr T fastInverseSqrt(T x, int iterations = 1) noexcept {
static_assert(std::numeric_limits<T>::is_iec559 &&
"Only IEEE 754 is supported");
static_assert(std::is_same_v<T, float> || std::is_same_v<T, double>,
"Only float and double are supported");
using I = std::conditional_t<std::is_same_v<T, float>, std::uint32_t,
std::uint64_t>;

constexpr I magic =
std::is_same_v<T, float> ? 0x5f3759df : 0x5fe6eb50c7b537a9;

union {
T f;
I i;
} u = {x};

u.i = magic - (u.i >> 1);

for (int i = 0; i < iterations; ++i) {
u.f *= 1.5f - 0.5f * x * u.f * u.f;
}

return u.f;
}

/// @brief Fast power function @see `std::pow`
/// Taken from
/// https://martin.ankerl.com/2012/01/25/optimized-approximative-pow-in-c-and-cpp
/// @param a the base
/// @param b the exponent
constexpr double fastPow(double a, double b) {
// enable only on IEEE 754
static_assert(std::numeric_limits<double>::is_iec559);

constexpr std::int64_t magic = 0x3FEF127F00000000;

union {
double f;
std::int64_t i;
} u = {a};

u.i = static_cast<std::int64_t>(b * (u.i - magic) + magic);

return u.f;
}

/// @brief Fast power function more precise than @see `fastPow`
/// Taken from
/// https://martin.ankerl.com/2012/01/25/optimized-approximative-pow-in-c-and-cpp
/// @param a the base
/// @param b the exponent
constexpr double fastPowMorePrecise(double a, double b) {
// binary exponentiation
double r = 1.0;
int exp = std::abs(static_cast<int>(b));
double base = b > 0 ? a : 1 / a;
while (exp != 0) {
if ((exp & 1) != 0) {
r *= base;
}
base *= base;
exp >>= 1;
}

return r * fastPow(a, b - static_cast<int>(b));
}

} // namespace Acts
6 changes: 5 additions & 1 deletion Examples/Python/src/Geometry.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -207,10 +207,14 @@ void addExperimentalGeometry(Context& ctx) {

// Detector volume definition
py::class_<DetectorVolume, std::shared_ptr<DetectorVolume>>(m,
"DetectorVolume");
"DetectorVolume")
.def("surfaces", &DetectorVolume::surfaces)
.def("surfacePtrs", &DetectorVolume::surfacePtrs);

// Detector definition
py::class_<Detector, std::shared_ptr<Detector>>(m, "Detector")
.def("volumes", &Detector::volumes)
.def("volumePtrs", &Detector::volumePtrs)
.def("numberVolumes",
[](Detector& self) { return self.volumes().size(); })
.def("extractMaterialSurfaces", [](Detector& self) {
Expand Down
17 changes: 14 additions & 3 deletions Examples/Scripts/Python/detector_creation.py
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,17 @@
geoContext = acts.GeometryContext()
[detector, contextors, store] = dd4hepDetector.finalize(geoContext, cOptions)

# OBJ style output
surfaces = []
for vol in detector.volumePtrs():
for surf in vol.surfacePtrs():
if surf.geometryId().sensitive() > 0:
surfaces.append(surf)
acts.examples.writeSurfacesObj(
surfaces, geoContext, [0, 120, 120], "odd-surfaces.obj"
)

# SVG style output
surfaceStyle = acts.svg.Style()
surfaceStyle.fillColor = [5, 150, 245]
surfaceStyle.fillOpacity = 0.5
Expand All @@ -58,7 +69,7 @@
volumeOptions = acts.svg.DetectorVolumeOptions()
volumeOptions.surfaceOptions = surfaceOptions

for ivol in range(detector.number_volumes()):
for ivol in range(detector.numberVolumes()):
acts.svg.viewDetector(
geoContext,
detector,
Expand All @@ -75,9 +86,9 @@
geoContext,
detector,
"odd",
[[ivol, volumeOptions] for ivol in range(detector.number_volumes())],
[[ivol, volumeOptions] for ivol in range(detector.numberVolumes())],
[["xy", ["sensitives"], xyRange], ["zr", ["materials"], zrRange]],
"detector",
)

# acts.examples.writeDetectorToJsonDetray(geoContext, detector, "odd-detray")
acts.examples.writeDetectorToJsonDetray(geoContext, detector, "odd-detray")
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,7 @@ void convert(nlohmann::json& jIndexedSurfaces,
jAccLink["type"] =
DetrayJsonHelper::accelerationLink(indexedSurfaces.casts);
jAccLink["index"] = std::numeric_limits<std::size_t>::max();
jIndexedSurfaces["acc_link"] = jAccLink;
jIndexedSurfaces["grid_link"] = jAccLink;
}
}
}
Expand Down
14 changes: 10 additions & 4 deletions Plugins/Json/src/DetectorJsonConverter.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -140,7 +140,7 @@ nlohmann::json Acts::DetectorJsonConverter::toJsonDetray(
// (2) surface grid section
nlohmann::json jSurfaceGrids;
nlohmann::json jSurfaceGridsData;
nlohmann::json jSurfaceGridsCollection;
nlohmann::json jSurfaceGridsInfoCollection;
nlohmann::json jSurfaceGridsHeader;
for (const auto [iv, volume] : enumerate(volumes)) {
// And its surface navigation delegates
Expand All @@ -150,15 +150,21 @@ nlohmann::json Acts::DetectorJsonConverter::toJsonDetray(
continue;
}
// Colplete the grid json for detray usage
jSurfacesDelegate["volume_link"] = iv;
jSurfacesDelegate["owner_link"] = iv;
// jSurfacesDelegate["acc_link"] =
nlohmann::json jSurfaceGridsCollection;
jSurfaceGridsCollection.push_back(jSurfacesDelegate);

nlohmann::json jSurfaceGridsInfo;
jSurfaceGridsInfo["volume_link"] = iv;
jSurfaceGridsInfo["grid_data"] = jSurfaceGridsCollection;
jSurfaceGridsInfoCollection.push_back(jSurfaceGridsInfo);
}
jSurfaceGridsData["grids"] = jSurfaceGridsCollection;
jSurfaceGridsData["grids"] = jSurfaceGridsInfoCollection;

jCommonHeader["tag"] = "surface_grids";
jSurfaceGridsHeader["common"] = jCommonHeader;
jSurfaceGridsHeader["grid_count"] = jSurfaceGridsCollection.size();
jSurfaceGridsHeader["grid_count"] = jSurfaceGridsInfoCollection.size();

jSurfaceGrids["header"] = jSurfaceGridsHeader;
jSurfaceGrids["data"] = jSurfaceGridsData;
Expand Down
1 change: 1 addition & 0 deletions Tests/Benchmarks/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -29,3 +29,4 @@ add_benchmark(SurfaceIntersection SurfaceIntersectionBenchmark.cpp)
add_benchmark(RayFrustum RayFrustumBenchmark.cpp)
add_benchmark(AnnulusBounds AnnulusBoundsBenchmark.cpp)
add_benchmark(StraightLineStepper StraightLineStepperBenchmark.cpp)
add_benchmark(QuickMath QuickMathBenchmark.cpp)
92 changes: 92 additions & 0 deletions Tests/Benchmarks/QuickMathBenchmark.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,92 @@
// This file is part of the Acts project.
//
// Copyright (C) 2024 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 <boost/test/data/test_case.hpp>
#include <boost/test/unit_test.hpp>

#include "Acts/Tests/CommonHelpers/BenchmarkTools.hpp"
#include "Acts/Utilities/QuickMath.hpp"

#include <cmath>
#include <random>

namespace bdata = boost::unit_test::data;

namespace Acts::Test {

/// @brief Another fast power function @see `fastPow`
// Taken from
// https://martin.ankerl.com/2007/02/11/optimized-exponential-functions-for-java
/// @param a the base
/// @param b the exponent
constexpr double fastPowAnother(double a, double b) {
// enable only on IEEE 754
static_assert(std::numeric_limits<double>::is_iec559);

union {
double f;
std::int64_t i;
} u = {};

u.i = static_cast<std::int64_t>(
9076650 * (a - 1) / (a + 1 + 4 * std::sqrt(a)) * b + 1072632447);
u.i <<= 32;

// result seems broken?
return u.f;
}

// Some randomness & number crunching
const unsigned int nTests = 10;
const unsigned int nReps = 10000;

BOOST_DATA_TEST_CASE(
benchmark_pow_25,
bdata::random(
(bdata::engine = std::mt19937(), bdata::seed = 21,
bdata::distribution = std::uniform_real_distribution<double>(-4, 4))) ^
bdata::xrange(nTests),
baseExp, index) {
(void)index;

const double base = std::pow(10, baseExp);
const double exp = 0.25;

std::cout << std::endl
<< "Benchmarking base=" << base << ", exp=" << exp << "..."
<< std::endl;
std::cout << "- void: "
<< Acts::Test::microBenchmark([&] { return 0; }, nReps)
<< std::endl;
std::cout << "- std::pow: "
<< Acts::Test::microBenchmark([&] { return std::pow(base, exp); },
nReps)
<< std::endl;
std::cout << "- std::exp: "
<< Acts::Test::microBenchmark(
[&] { return std::exp(std::log(base) * exp); }, nReps)
<< std::endl;
std::cout << "- std::sqrt: "
<< Acts::Test::microBenchmark(
[&] { return std::sqrt(std::sqrt(base)); }, nReps)
<< std::endl;
std::cout << "- fastPow: "
<< Acts::Test::microBenchmark([&] { return fastPow(base, exp); },
nReps)
<< std::endl;
std::cout << "- fastPowMorePrecise: "
<< Acts::Test::microBenchmark(
[&] { return fastPowMorePrecise(base, exp); }, nReps)
<< std::endl;
std::cout << "- fastPowAnother: "
<< Acts::Test::microBenchmark(
[&] { return fastPowAnother(base, exp); }, nReps)
<< std::endl;
}

} // namespace Acts::Test
1 change: 1 addition & 0 deletions Tests/UnitTests/Core/Utilities/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -51,3 +51,4 @@ target_compile_definitions(ActsUnitTestAnyDebug PRIVATE _ACTS_ANY_ENABLE_VERBOSE
add_unittest(ParticleData ParticleDataTests.cpp)
add_unittest(Zip ZipTests.cpp)
add_unittest(TransformRange TransformRangeTests.cpp)
add_unittest(QuickMath QuickMathTests.cpp)
Loading

0 comments on commit 3e6442c

Please sign in to comment.