From b10bbc6171b88c799f4e4b31afa6365730c52976 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 19 May 2026 17:13:17 -0400 Subject: [PATCH 01/28] Allow zeros and optimized poisson distribution for each constructed case --- .../distribution/PoissonDistribution.hh | 83 +++++++++++++------ .../distribution/PoissonDistribution.test.cc | 19 +++++ 2 files changed, 77 insertions(+), 25 deletions(-) diff --git a/src/corecel/random/distribution/PoissonDistribution.hh b/src/corecel/random/distribution/PoissonDistribution.hh index 0d0639d124..68ac6949f0 100644 --- a/src/corecel/random/distribution/PoissonDistribution.hh +++ b/src/corecel/random/distribution/PoissonDistribution.hh @@ -9,10 +9,8 @@ #include #include "corecel/Assert.hh" -#include "corecel/Constants.hh" #include "corecel/Macros.hh" #include "corecel/Types.hh" -#include "corecel/math/Algorithms.hh" #include "NormalDistribution.hh" @@ -20,7 +18,7 @@ namespace celeritas { //---------------------------------------------------------------------------// /*! - * Sample from a Poisson distribution. + * Sample from a generalized Poisson distribution. * * The Poisson distribution describes the probability of \f$ k \f$ events * occurring in a fixed interval given a mean rate of occurrence \f$ \lambda @@ -45,9 +43,13 @@ namespace celeritas * approximation for \f$ \lambda > 16 \f$ (see \c G4Poisson), which is faster * but less accurate than other methods. The same approach is used here. * - * \todo Break this into two distributions: one actual poisson distribution, - * one "integer normal" distribution, and a variant type that selects between - * them. In most cases we care about, lambda is small. + * In the degenerate case of \f$ \lambda = 0 \f$, the result is always zero. + * + * \todo Rename to GeneralPoissonDistribution (or something similar) since + * it's effectively an inefficient variant combining: + * - an actual poisson distribution, + * - an "integer normal" distribution, and + * - a "zero" distribution. */ template class PoissonDistribution @@ -60,8 +62,11 @@ class PoissonDistribution //!@} public: - // Construct with defaults - explicit inline CELER_FUNCTION PoissonDistribution(real_type lambda = 1); + // Construct with distribution parameter + explicit inline CELER_FUNCTION PoissonDistribution(real_type lambda); + + //! Construct with default lambda of 1 + CELER_FUNCTION PoissonDistribution() : PoissonDistribution{1} {} // Sample a random number according to the distribution template @@ -71,8 +76,15 @@ class PoissonDistribution static CELER_CONSTEXPR_FUNCTION int lambda_threshold() { return 16; } private: - real_type const lambda_; - NormalDistribution sample_normal_; + enum class Method + { + zero, + poisson, + gaussian + }; + Method method_; + real_type exp_lambda_{}; + NormalDistribution sample_normal_{}; }; //---------------------------------------------------------------------------// @@ -84,9 +96,24 @@ class PoissonDistribution template CELER_FUNCTION PoissonDistribution::PoissonDistribution(real_type lambda) - : lambda_(lambda), sample_normal_(lambda_, std::sqrt(lambda_)) { - CELER_EXPECT(lambda_ > 0); + CELER_EXPECT(lambda >= 0); + + if (lambda <= 0) + { + method_ = Method::zero; + } + else if (lambda <= PoissonDistribution::lambda_threshold()) + { + method_ = Method::poisson; + exp_lambda_ = std::exp(lambda); + } + else + { + method_ = Method::gaussian; + // Add 0.5 to mean for correct rounding + sample_normal_ = NormalDistribution{lambda + 0.5, std::sqrt(lambda)}; + } } //---------------------------------------------------------------------------// @@ -98,20 +125,26 @@ template CELER_FUNCTION auto PoissonDistribution::operator()(Generator& rng) -> result_type { - if (lambda_ <= PoissonDistribution::lambda_threshold()) + switch (method_) { - // Use direct method - int k = 0; - real_type p = std::exp(lambda_); - do - { - ++k; - p *= generate_canonical(rng); - } while (p > 1); - return static_cast(k - 1); - } - // Use Gaussian approximation rounded to nearest integer - return result_type(sample_normal_(rng) + real_type(0.5)); + case Method::zero: + return 0; + case Method::poisson: { + int k = 0; + real_type p = exp_lambda_; + do + { + ++k; + p *= generate_canonical(rng); + } while (p > 1); + return static_cast(k - 1); + } + case Method::gaussian: + // Use Gaussian approximation rounded to nearest integer + return static_cast( + clamp_to_nonneg(sample_normal_(rng))); + }; } + //---------------------------------------------------------------------------// } // namespace celeritas diff --git a/test/corecel/random/distribution/PoissonDistribution.test.cc b/test/corecel/random/distribution/PoissonDistribution.test.cc index dde60a2339..b76b224b1b 100644 --- a/test/corecel/random/distribution/PoissonDistribution.test.cc +++ b/test/corecel/random/distribution/PoissonDistribution.test.cc @@ -18,6 +18,25 @@ namespace test { //---------------------------------------------------------------------------// +TEST(PoissonDistributionTest, bin_zero) +{ + int num_samples = 100; + + // Small lambda will use the direct method, which requires on average + // lambda + 1 RNG samples + PoissonDistribution sample_poisson{0}; + DiagnosticRngEngine rng; + + Histogram histogram(4, {0, 1e-3}); + for ([[maybe_unused]] int i : range(num_samples)) + { + histogram(sample_poisson(rng)); + } + static unsigned int const expected_counts[] = {100, 0, 0, 0}; + EXPECT_VEC_EQ(expected_counts, histogram.counts()); + EXPECT_EQ(0, rng.count()); +} + TEST(PoissonDistributionTest, bin_small) { int num_samples = 10000; From f4748f0c662d2f144d3aab0ef0d9fea58a22e0d5 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 19 May 2026 17:14:50 -0400 Subject: [PATCH 02/28] Remove now-unnecessary edge cases --- .../distribution/EnergyLossUrbanDistribution.hh | 2 +- src/celeritas/optical/gen/CherenkovOffload.hh | 17 ++++++----------- 2 files changed, 7 insertions(+), 12 deletions(-) diff --git a/src/celeritas/em/distribution/EnergyLossUrbanDistribution.hh b/src/celeritas/em/distribution/EnergyLossUrbanDistribution.hh index 0bfac1eaf1..29527c2db9 100644 --- a/src/celeritas/em/distribution/EnergyLossUrbanDistribution.hh +++ b/src/celeritas/em/distribution/EnergyLossUrbanDistribution.hh @@ -270,7 +270,7 @@ EnergyLossUrbanDistribution::sample_excitation_loss(Engine& rng) mean += xs_exc_[i] * binding_energy_[i]; variance += xs_exc_[i] * ipow<2>(binding_energy_[i]); } - else if (xs_exc_[i] > 0) + else { // The loss due to excitation is \f$ \Delta E_{exc} = n_1 E_1 + n_2 // E_2 \f$, where the number of collisions \f$ n_i \f$ is sampled diff --git a/src/celeritas/optical/gen/CherenkovOffload.hh b/src/celeritas/optical/gen/CherenkovOffload.hh index 96f01741c8..f252b7ceca 100644 --- a/src/celeritas/optical/gen/CherenkovOffload.hh +++ b/src/celeritas/optical/gen/CherenkovOffload.hh @@ -67,7 +67,7 @@ class CherenkovOffload real_type step_length_; OffloadPreStepData const& pre_step_; optical::GeneratorStepData post_step_; - real_type num_photons_per_len_; + PoissonDistribution sample_num_photons_{0}; }; //---------------------------------------------------------------------------// @@ -96,11 +96,12 @@ CherenkovOffload::CherenkovOffload(units::ElementaryCharge charge, using namespace celeritas::literals; - units::LightSpeed beta( + // NOTE: this uses the average beta, not the average number of photons + units::LightSpeed avg_beta( 0.5_r * (pre_step_.speed.value() + post_step_.speed.value())); - CherenkovDndxCalculator calculate_dndx(pre_mat, shared, charge_); - num_photons_per_len_ = calculate_dndx(beta); + sample_num_photons_ + = PoissonDistribution{calculate_dndx(avg_beta) * step_length_}; } //---------------------------------------------------------------------------// @@ -135,14 +136,8 @@ template CELER_FUNCTION optical::GeneratorDistributionData CherenkovOffload::operator()(Generator& rng) { - if (num_photons_per_len_ == 0) - { - return {}; - } - optical::GeneratorDistributionData data; - data.num_photons = PoissonDistribution(num_photons_per_len_ - * step_length_)(rng); + data.num_photons = sample_num_photons_(rng); if (data.num_photons > 0) { data.type = GeneratorType::cherenkov; From b05f48995afda4ed425f6ebbde4b4d60b2956cfc Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 19 May 2026 17:15:30 -0400 Subject: [PATCH 03/28] Slightly optimize gaussian-possion-like calculation --- src/celeritas/optical/gen/ScintillationOffload.hh | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/celeritas/optical/gen/ScintillationOffload.hh b/src/celeritas/optical/gen/ScintillationOffload.hh index 5aae2355b2..46f39e3426 100644 --- a/src/celeritas/optical/gen/ScintillationOffload.hh +++ b/src/celeritas/optical/gen/ScintillationOffload.hh @@ -131,9 +131,9 @@ ScintillationOffload::operator()(Generator& rng) real_type sigma = shared_.resolution_scale[pre_step_.material] * std::sqrt(mean_num_photons_); - result.num_photons = static_cast(clamp_to_nonneg( - NormalDistribution(mean_num_photons_, sigma)(rng) - + 0.5_r)); + result.num_photons = static_cast( + clamp_to_nonneg(NormalDistribution( + mean_num_photons_ + 0.5_r, sigma)(rng))); } else if (mean_num_photons_ > 0) { From 624bb5123972c7060847a633cdd94bffdae4f682 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Thu, 21 May 2026 06:51:59 -0400 Subject: [PATCH 04/28] Add notes and zero-sample poisson for WLS interactor --- .../interactor/WavelengthShiftInteractor.hh | 45 ++++++++++--------- 1 file changed, 24 insertions(+), 21 deletions(-) diff --git a/src/celeritas/optical/interactor/WavelengthShiftInteractor.hh b/src/celeritas/optical/interactor/WavelengthShiftInteractor.hh index 86439192d7..ccf6e05412 100644 --- a/src/celeritas/optical/interactor/WavelengthShiftInteractor.hh +++ b/src/celeritas/optical/interactor/WavelengthShiftInteractor.hh @@ -10,12 +10,10 @@ #include "corecel/Types.hh" #include "corecel/random/distribution/PoissonDistribution.hh" #include "celeritas/Types.hh" -#include "celeritas/grid/NonuniformGridCalculator.hh" #include "celeritas/optical/Interaction.hh" #include "celeritas/optical/ParticleTrackView.hh" #include "celeritas/optical/SimTrackView.hh" #include "celeritas/optical/WavelengthShiftData.hh" -#include "celeritas/phys/InteractionUtils.hh" namespace celeritas { @@ -30,6 +28,11 @@ namespace optical * * \todo See if initializing the first photon directly in this track slot * improves performance + * \todo This is shoehorned into the \c Interactor paradigm: replace with + * \c WlsOffloadSampler that samples a new \c WlsDistribution . (If we replace + * the multiple WLS models with a single WLS model, similar to how + * scintillation has multiple time and energy spectra, this could also sample + * the spectrum ID.) */ class WavelengthShiftInteractor { @@ -87,12 +90,22 @@ WavelengthShiftInteractor::WavelengthShiftInteractor( CELER_EXPECT(distribution_id_ < data_.distributions.size()); CELER_EXPECT(!data_.distributions[distribution_id_]); - distribution_.type = shared.type; - distribution_.energy = particle.energy(); - distribution_.time = sim.time(); - distribution_.position = pos; - distribution_.primary = sim.primary_id(); - distribution_.material = mat_id; + if (distribution_.energy <= emission_threshold_) + { + // If the incident particle energy is below the lower bound of the + // emitted energy sampling grid, don't emit any photons + sample_num_photons_ = PoissonDistribution<>{0}; + distribution_.num_photons = 0; + } + else + { + distribution_.type = shared.type; + distribution_.energy = particle.energy(); + distribution_.time = sim.time(); + distribution_.position = pos; + distribution_.primary = sim.primary_id(); + distribution_.material = mat_id; + } } //---------------------------------------------------------------------------// @@ -106,20 +119,10 @@ CELER_FUNCTION Interaction WavelengthShiftInteractor::operator()(Engine& rng) { Interaction result = Interaction::from_absorption(); - if (distribution_.energy <= emission_threshold_) - { - // If the incident particle energy is below the lower bound of the - // emitted energy sampling grid, don't emit any photons - distribution_.num_photons = 0; - } - else - { - // Sample the number of photons generated from WLS. - distribution_.num_photons = sample_num_photons_(rng); - } - CELER_ASSERT(distribution_ || distribution_.num_photons == 0); + // Sample the number of photons generated from WLS + distribution_.num_photons = sample_num_photons_(rng); data_.distributions[distribution_id_] = distribution_; - + CELER_ENSURE(distribution_ || distribution_.num_photons == 0); return result; } From 3b53868c59222ba57cbb9b7e1b18d0149deacc2a Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Thu, 21 May 2026 11:46:32 -0400 Subject: [PATCH 05/28] Address feedback --- src/celeritas/optical/interactor/WavelengthShiftInteractor.hh | 2 +- src/corecel/random/distribution/PoissonDistribution.hh | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/celeritas/optical/interactor/WavelengthShiftInteractor.hh b/src/celeritas/optical/interactor/WavelengthShiftInteractor.hh index ccf6e05412..60c22d4b40 100644 --- a/src/celeritas/optical/interactor/WavelengthShiftInteractor.hh +++ b/src/celeritas/optical/interactor/WavelengthShiftInteractor.hh @@ -90,7 +90,7 @@ WavelengthShiftInteractor::WavelengthShiftInteractor( CELER_EXPECT(distribution_id_ < data_.distributions.size()); CELER_EXPECT(!data_.distributions[distribution_id_]); - if (distribution_.energy <= emission_threshold_) + if (particle.energy() <= emission_threshold_) { // If the incident particle energy is below the lower bound of the // emitted energy sampling grid, don't emit any photons diff --git a/src/corecel/random/distribution/PoissonDistribution.hh b/src/corecel/random/distribution/PoissonDistribution.hh index 68ac6949f0..168f381423 100644 --- a/src/corecel/random/distribution/PoissonDistribution.hh +++ b/src/corecel/random/distribution/PoissonDistribution.hh @@ -11,6 +11,7 @@ #include "corecel/Assert.hh" #include "corecel/Macros.hh" #include "corecel/Types.hh" +#include "corecel/math/Algorithms.hh" #include "NormalDistribution.hh" From 9886120ac9f8f4fa65e28253adc9d53c86c9fe4c Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Fri, 22 May 2026 07:37:27 -0400 Subject: [PATCH 06/28] Add distribution record (but not generic support) for uniform real to allow testing --- src/corecel/random/data/DistributionData.hh | 11 ++++++++++- .../random/distribution/UniformRealDistribution.hh | 3 ++- 2 files changed, 12 insertions(+), 2 deletions(-) diff --git a/src/corecel/random/data/DistributionData.hh b/src/corecel/random/data/DistributionData.hh index efbddb2f0f..b23186d930 100644 --- a/src/corecel/random/data/DistributionData.hh +++ b/src/corecel/random/data/DistributionData.hh @@ -7,7 +7,6 @@ #pragma once #include "corecel/Macros.hh" -#include "corecel/OpaqueId.hh" #include "corecel/Types.hh" #include "corecel/cont/Array.hh" #include "corecel/data/Collection.hh" @@ -38,6 +37,16 @@ struct DeltaDistributionRecord T value{}; }; +//---------------------------------------------------------------------------// +/*! + * Data for sampling a value from a uniform distribution. + */ +struct UniformRealDistributionRecord +{ + real_type lower{}; + real_type upper{}; +}; + //---------------------------------------------------------------------------// /*! * Data for sampling from a normal distribution. diff --git a/src/corecel/random/distribution/UniformRealDistribution.hh b/src/corecel/random/distribution/UniformRealDistribution.hh index 3f2195678e..5745a39236 100644 --- a/src/corecel/random/distribution/UniformRealDistribution.hh +++ b/src/corecel/random/distribution/UniformRealDistribution.hh @@ -9,9 +9,9 @@ #include #include -#include "corecel/Assert.hh" #include "corecel/Macros.hh" #include "corecel/Types.hh" +#include "corecel/random/data/DistributionData.hh" #include "GenerateCanonical.hh" @@ -43,6 +43,7 @@ class UniformRealDistribution //! \name Type aliases using real_type = RealType; using result_type = real_type; + using RecordT = UniformRealDistributionRecord; //!@} public: From 7b2d0e6e8e0c233680e3623786e839efae30ddc8 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Fri, 22 May 2026 08:50:57 -0400 Subject: [PATCH 07/28] Add rejection assertion checking to TruncatedDistribution --- .../distribution/TruncatedDistribution.hh | 20 +++++++++++- .../TruncatedDistribution.test.cc | 31 +++++++++++++++++++ 2 files changed, 50 insertions(+), 1 deletion(-) diff --git a/src/corecel/random/distribution/TruncatedDistribution.hh b/src/corecel/random/distribution/TruncatedDistribution.hh index 79d047e610..59244cf788 100644 --- a/src/corecel/random/distribution/TruncatedDistribution.hh +++ b/src/corecel/random/distribution/TruncatedDistribution.hh @@ -6,9 +6,10 @@ //---------------------------------------------------------------------------// #pragma once +#include "corecel/Config.hh" + #include "corecel/Assert.hh" #include "corecel/Macros.hh" -#include "corecel/Types.hh" #include "corecel/random/data/DistributionData.hh" namespace celeritas @@ -20,6 +21,9 @@ namespace celeritas * Sample from an arbitrary distribution truncated to a finite interval. Values * are drawn from the underlying distribution using rejection sampling until * they fall within the truncation bounds. + * + * \warning Because of the rejection sampling it is possible to create a + * distribution */ template class TruncatedDistribution @@ -49,6 +53,9 @@ class TruncatedDistribution template inline CELER_FUNCTION result_type operator()(Generator& rng); + // Assert fewer than this number of samples is tried (CELERITAS_DEBUG only) + static constexpr inline int max_debug_samples = 32; + private: Distribution sample_; real_type lower_; @@ -81,11 +88,22 @@ template CELER_FUNCTION auto TruncatedDistribution::operator()(Generator& rng) -> result_type { + int num_remaining_samples = max_debug_samples; result_type result; do { // Reject samples outside the truncation bounds result = sample_(rng); + // Prevent infinite loops (debug assertions only) + if constexpr (CELERITAS_DEBUG) + { + if (--num_remaining_samples < 0) + { + CELER_DEBUG_FAIL( + "too many samples taken in TruncatedDistribution", + internal); + } + } } while (result < lower_ || result > upper_); return result; } diff --git a/test/corecel/random/distribution/TruncatedDistribution.test.cc b/test/corecel/random/distribution/TruncatedDistribution.test.cc index fa481d9591..80231c0bfd 100644 --- a/test/corecel/random/distribution/TruncatedDistribution.test.cc +++ b/test/corecel/random/distribution/TruncatedDistribution.test.cc @@ -8,10 +8,13 @@ #include +#include "corecel/Assert.hh" #include "corecel/cont/Range.hh" #include "corecel/random/DiagnosticRngEngine.hh" #include "corecel/random/Histogram.hh" +#include "corecel/random/HistogramSampler.hh" #include "corecel/random/distribution/NormalDistribution.hh" +#include "corecel/random/distribution/UniformRealDistribution.hh" #include "celeritas_test.hh" @@ -47,6 +50,34 @@ TEST(TruncatedDistributionTest, normal) EXPECT_EQ(20984, rng.count()); } +TEST(TruncatedDistributionTest, uniform) +{ + using TruncatedUniformDbl + = TruncatedDistribution>; + + // Check that all values are within truncated range + HistogramSampler calc_histogram(8, {0, 4}, 100); + // Underlying distribution samples between -1, 4 + auto result = calc_histogram(TruncatedUniformDbl(0, 4, -1, 4)); + + SampledHistogram ref; + ref.distribution = {0.16, 0.26, 0.28, 0.26, 0.14, 0.34, 0.24, 0.32}; + ref.rng_count = 2.36; + EXPECT_REF_EQ(ref, result); +} + +TEST(TruncatedDistributionTest, TEST_IF_CELERITAS_DEBUG(rejection)) +{ + using TruncatedUniformDbl + = TruncatedDistribution>; + // Truncated bounds are outside actual distribution + TruncatedUniformDbl sample_forever(-10, -1, 1, 5); + + DiagnosticRngEngine rng; + EXPECT_THROW(sample_forever(rng), DebugError); + EXPECT_EQ(66, rng.exchange_count()); +} + //---------------------------------------------------------------------------// } // namespace test } // namespace celeritas From 24943789537eab618f8cccc5d29f59eb0520a9fc Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Fri, 22 May 2026 09:08:23 -0400 Subject: [PATCH 08/28] Add documentation to testing distributions --- doc/CMakeLists.txt | 2 + doc/development/testing.rst | 64 +++++++++++++------ .../distribution/TruncatedDistribution.hh | 12 +++- test/Test.hh | 16 ++++- test/corecel/ScopedLogStorer.hh | 1 + test/corecel/StringSimplifier.hh | 14 +++- test/corecel/random/DiagnosticRngEngine.hh | 2 +- test/corecel/random/Histogram.hh | 1 + test/corecel/random/HistogramSampler.hh | 12 +++- 9 files changed, 92 insertions(+), 32 deletions(-) diff --git a/doc/CMakeLists.txt b/doc/CMakeLists.txt index a238181a62..c0edef34a8 100644 --- a/doc/CMakeLists.txt +++ b/doc/CMakeLists.txt @@ -288,6 +288,8 @@ file(GLOB _DOXYGEN_SOURCE "${PROJECT_SOURCE_DIR}/src/celeritas/optical/surface/*.hh" "${PROJECT_SOURCE_DIR}/src/celeritas/optical/surface/model/*.hh" "${PROJECT_SOURCE_DIR}/test/*.hh" + "${PROJECT_SOURCE_DIR}/test/corecel/*.hh" + "${PROJECT_SOURCE_DIR}/test/corecel/random/*.hh" ) # IMPORTANT: the target name and USE_STAMP_FILE generate a dependency used by diff --git a/doc/development/testing.rst b/doc/development/testing.rst index 4041c21036..e0048c0014 100644 --- a/doc/development/testing.rst +++ b/doc/development/testing.rst @@ -20,32 +20,21 @@ functions to better isolate conditionals from each other. .. _cyclomatic complexity: https://en.wikipedia.org/wiki/Cyclomatic_complexity -Running CTest -------------- - -When configured with ``CELERITAS_BUILD_TESTS`` (see :ref:`configuration`), -CTest_ will be automatically configured. Running through CTest sets special -environment variables for data, disabling GPUs, or testing the code itself. -CTest can be run either through ``ninja test`` or by manually invoking -``ctest``. Two useful ways to run are ``ctest -V -R ``, which will -run one or more tests that match the given regular expression (such as -``corecel/math/``), and ``ctest -j --output-on-failure`` which runs in parallel -and prints only test failures. - -.. _CTest: https://cmake.org/cmake/help/latest/manual/ctest.1.html - -Using GoogleTest ----------------- +Writing unit tests +------------------ `Google test`_ is very well documented, and because so much testing code exists in AI training data, tools like ChatGPT and Copilot are very good at writing a first pass at test code. + Celeritas defines a base class test harness with some utility functions: .. doxygenclass:: celeritas::test::Test -as well as several macros that simplify testing with floating-point data (and -vectors thereof): +Testing real numbers and JSON data +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Celeritas defines special test macros that simplify testing with floating-point data (and vectors thereof): .. doxygendefine:: EXPECT_VEC_EQ .. doxygendefine:: EXPECT_REAL_EQ @@ -53,15 +42,50 @@ vectors thereof): .. doxygendefine:: EXPECT_SOFT_NEAR .. doxygendefine:: EXPECT_VEC_SOFT_EQ .. doxygendefine:: EXPECT_VEC_NEAR +.. doxygendefine:: EXPECT_REF_EQ .. doxygendefine:: EXPECT_JSON_EQ + +Expected data can be printed with another macro: + .. doxygendefine:: PRINT_EXPECTED For more details on the test harnesses, especially the hierarchy used for setting up physics problems for testing, see the ``celeritas::test`` namespace in the Doxygen developer documentation. -You can run most tests manually from the build directory and filter so that -only a subset of tests run:: +Testing random number distributions +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +.. doxygenclass:: celeritas::test::DiagnosticRngEngine + +.. doxygenclass:: celeritas::test::Histogram + +.. doxygenclass:: celeritas::test::HistogramSampler + +Testing strings and logger output +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +.. doxygenclass:: celeritas::test::ScopedLogStorer + +.. doxygenclass:: celeritas::test::StringSimplifier + + +Running unit tests +------------------ + +When configured with ``CELERITAS_BUILD_TESTS`` (see :ref:`configuration`), +CTest_ will be automatically configured. Running through CTest sets special +environment variables for data, disabling GPUs, or testing the code itself. +CTest can be run either through ``ninja test`` or by manually invoking +``ctest``. Two useful ways to run are ``ctest -V -R ``, which will +run one or more tests that match the given regular expression (such as +``corecel/math/``), and ``ctest -j --output-on-failure`` which runs in parallel +and prints only test failures. + +.. _CTest: https://cmake.org/cmake/help/latest/manual/ctest.1.html + +You can *also* run most tests manually from the build directory and filter so +that only a subset of tests run:: $ ./test/celeritas/global_Stepper --gtest_filter=SimpleComptonTest.host diff --git a/src/corecel/random/distribution/TruncatedDistribution.hh b/src/corecel/random/distribution/TruncatedDistribution.hh index 59244cf788..273dbc95ac 100644 --- a/src/corecel/random/distribution/TruncatedDistribution.hh +++ b/src/corecel/random/distribution/TruncatedDistribution.hh @@ -22,8 +22,14 @@ namespace celeritas * are drawn from the underlying distribution using rejection sampling until * they fall within the truncation bounds. * - * \warning Because of the rejection sampling it is possible to create a - * distribution + * \warning Because of the rejection sampling, it is possible to create a + * distribution that never (or almost never) accepts a value from the + * underlying distribution. A loop checker is active in debug mode that + * prevents more than 32 failed samples. Since GPU performance is so sensitive + * to rejection failures, consider transforming the underlying distribution if + * this limit is hit. When using this distribution for physics distributions, + * it is is \em strongly recommended to use \c test::HistogramSampler to verify + * the number of samples being taken. */ template class TruncatedDistribution @@ -53,7 +59,7 @@ class TruncatedDistribution template inline CELER_FUNCTION result_type operator()(Generator& rng); - // Assert fewer than this number of samples is tried (CELERITAS_DEBUG only) + //! Assert fewer than this number of samples is tried when CELERITAS_DEBUG static constexpr inline int max_debug_samples = 32; private: diff --git a/test/Test.hh b/test/Test.hh index 0e8210d74e..25f352c514 100644 --- a/test/Test.hh +++ b/test/Test.hh @@ -23,6 +23,10 @@ namespace celeritas * By being a sub-namspace, no \c using declarations are needed to access * Celeritas functionality, and no symbol conflicts with main Celeritas code * will accidentally occur. + * + * Tests of the \c detail (or other) namespace functionality should be nested + * in a separate namespace \c celeritas::detail::test . They will have to + * explicitly use \c celeritas::test::Test to access the \c Test class. */ namespace test { @@ -32,16 +36,24 @@ namespace test * * The test harness is constructed and destroyed once per subtest. It contains * helper functions and data commonly needed in Celeritas tests. + * + * - \c make_unique_filename : construct a filename unique to the given test + * suite and function, which will prevent files from stepping on each other + * when running tests in parallel + * - \c test_data_path : get the path to a test file at + * `{source}/test/{subdir}/data/{filename}` . It verifies that the path is an + * existing file before passing it to the code to be tested. + * - \c inf , \c inff: double- and single-precision infinities */ class Test : public ::testing::Test { public: Test() = default; - // Generate test-unique filename + //! Generate test-unique filename virtual std::string make_unique_filename(std::string_view ext); - // Make a unique filename with no extension + //! Make a unique filename with no extension std::string make_unique_filename() { return this->make_unique_filename({}); diff --git a/test/corecel/ScopedLogStorer.hh b/test/corecel/ScopedLogStorer.hh index 617329f64a..01e439f638 100644 --- a/test/corecel/ScopedLogStorer.hh +++ b/test/corecel/ScopedLogStorer.hh @@ -30,6 +30,7 @@ namespace test * You can use the \c CELER_LOG_SCOPED environment variable to print the * captured log messages as they are written. * + * \par Example: * \code ScopedLogStorer scoped_log_{&celeritas::world_logger()}; CELER_LOG(info) << "captured"; diff --git a/test/corecel/StringSimplifier.hh b/test/corecel/StringSimplifier.hh index 41fa812317..412409676c 100644 --- a/test/corecel/StringSimplifier.hh +++ b/test/corecel/StringSimplifier.hh @@ -21,19 +21,27 @@ namespace test * - Removes pointers * - Removes ANSI escape sequences * - Rounds floating points to a given digit of precision + * + * \par Example: + * \code + * auto simplify_str = StringSimplifier{3}; + * EXPECT_EQ("1.23 [s]", simplify_str(exception_msg)); + * \endcode */ class StringSimplifier { public: // Construct with number of digit of precision - inline StringSimplifier(int precision); - StringSimplifier() = default; + explicit inline StringSimplifier(int precision); + + //! Default to 4 digits of precision + StringSimplifier() : StringSimplifier{4} {} // Simplify [[nodiscard]] std::string operator()(std::string const& old) const; private: - int precision_{4}; + int precision_{}; std::string simplify_sci(std::string s) const; std::string simplify_float(std::string s) const; diff --git a/test/corecel/random/DiagnosticRngEngine.hh b/test/corecel/random/DiagnosticRngEngine.hh index 95c5ccc7d2..fa93394794 100644 --- a/test/corecel/random/DiagnosticRngEngine.hh +++ b/test/corecel/random/DiagnosticRngEngine.hh @@ -15,7 +15,7 @@ namespace test { //---------------------------------------------------------------------------// /*! - * Diagnostic wrapper that counts the number of calls to operator(). + * RNG wrapper that counts the number of calls to operator(). * * This wraps a low-level pseudorandom generator's call function. It can be * used to determine the efficiency of rejection algorithms without changing diff --git a/test/corecel/random/Histogram.hh b/test/corecel/random/Histogram.hh index dd8b76c953..523f970494 100644 --- a/test/corecel/random/Histogram.hh +++ b/test/corecel/random/Histogram.hh @@ -28,6 +28,7 @@ namespace test * overflow bins. All bins are half-open except for the rightmost bin, which * will include values equal to the upper domain boundary. * + * \par Example: * To test that all samples are within the domain: * \code EXPECT_EQ(0, hist.underflow()) diff --git a/test/corecel/random/HistogramSampler.hh b/test/corecel/random/HistogramSampler.hh index c710d54b88..e07b0060a2 100644 --- a/test/corecel/random/HistogramSampler.hh +++ b/test/corecel/random/HistogramSampler.hh @@ -29,19 +29,27 @@ struct SampledHistogram //! Average number of RNG samples double rng_count{}; + // Print expected code to cout void print_expected() const; + + // Print to a stream + friend std::ostream& operator<<(std::ostream&, SampledHistogram const&); }; //---------------------------------------------------------------------------// /*! * Sample one or more distributions, returning a histogram. * + * The sampled histogram is a \em density , so it is recommended to make the + * sampled width (delta of second parameter) \em and the number of samples + * (third parameter) evenly divisible into a power of 10 for prettier printing. + * + * \par Example: * \code constexpr size_type num_samples = 1000; HistogramSampler calc_histogram(8, {-1, 1}, num_samples); std::vector actual; - for (real_type inc_e : {0.1, 1.0, 1e2, 1e3, 1e6}) { for (real_type eps : {0.001, 0.01, 0.1}) @@ -91,8 +99,6 @@ class HistogramSampler // FREE FUNCTIONS //---------------------------------------------------------------------------// -std::ostream& operator<<(std::ostream& os, SampledHistogram const& sh); - ::testing::AssertionResult IsRefEq(char const* expr1, char const* expr2, SampledHistogram const& val1, From 55d9d2421c53d0330c3e9ea63eaafc0d7ef5dd4e Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Fri, 22 May 2026 09:26:20 -0400 Subject: [PATCH 09/28] Allow truncated distribution to be constructed via constants --- .../optical/surface/GaussianRoughnessSampler.hh | 2 +- .../random/distribution/TruncatedDistribution.hh | 10 +++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/celeritas/optical/surface/GaussianRoughnessSampler.hh b/src/celeritas/optical/surface/GaussianRoughnessSampler.hh index c56e5636ad..ff9e83db1a 100644 --- a/src/celeritas/optical/surface/GaussianRoughnessSampler.hh +++ b/src/celeritas/optical/surface/GaussianRoughnessSampler.hh @@ -74,7 +74,7 @@ CELER_FUNCTION GaussianRoughnessSampler::GaussianRoughnessSampler(Real3 const& normal, real_type sigma_alpha) : normal_(normal) - , sample_alpha_(-0.5 * constants::pi, 0.5 * constants::pi, 0, sigma_alpha) + , sample_alpha_(-constants::pi / 2, constants::pi / 2, 0, sigma_alpha) , f_max_(fmin(real_type{1}, 4 * sigma_alpha)) , reject_(f_max_) { diff --git a/src/corecel/random/distribution/TruncatedDistribution.hh b/src/corecel/random/distribution/TruncatedDistribution.hh index 273dbc95ac..40a8cb7fa3 100644 --- a/src/corecel/random/distribution/TruncatedDistribution.hh +++ b/src/corecel/random/distribution/TruncatedDistribution.hh @@ -44,9 +44,9 @@ class TruncatedDistribution public: // Construct with distribution and truncation bounds - template + template inline CELER_FUNCTION - TruncatedDistribution(real_type lower, real_type upper, Args&&... args); + TruncatedDistribution(T lower, U upper, Args&&... args); // Construct from record explicit CELER_FUNCTION TruncatedDistribution(RecordT const& record) @@ -75,10 +75,10 @@ class TruncatedDistribution * Construct with distribution and truncation bounds. */ template -template +template CELER_FUNCTION -TruncatedDistribution::TruncatedDistribution(real_type lower, - real_type upper, +TruncatedDistribution::TruncatedDistribution(T lower, + U upper, Args&&... args) : sample_(celeritas::forward(args)...), lower_(lower), upper_(upper) { From cb3f4c4cce5265768019d98a4e6edd5ecaab865d Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Fri, 22 May 2026 09:32:32 -0400 Subject: [PATCH 10/28] Add rounded nonnegative distribution wrapper Implement a templated RoundedNonnegDistribution that wraps an underlying floating-point distribution and returns rounded non-negative integers. Add focused unit tests for rounding/clamping behavior, one-underlying-sample usage, and deterministic uniform sampling output. Register the new test in corecel CMake test definitions. Prompt: "Using #sym:TruncatedDistribution as a template, create a RoundedNonnegDistribution that efficiently samples from an underlying distribution and returns a non-negative integer. Template parameters should be underlying distribution and integer type (default to celeritas::size_type)." Assisted-by: GitHub Copilot (GPT-5.3-Codex) --- .../distribution/RoundedNonnegDistribution.hh | 88 ++++++++++++ test/corecel/CMakeLists.txt | 1 + .../RoundedNonnegDistribution.test.cc | 132 ++++++++++++++++++ 3 files changed, 221 insertions(+) create mode 100644 src/corecel/random/distribution/RoundedNonnegDistribution.hh create mode 100644 test/corecel/random/distribution/RoundedNonnegDistribution.test.cc diff --git a/src/corecel/random/distribution/RoundedNonnegDistribution.hh b/src/corecel/random/distribution/RoundedNonnegDistribution.hh new file mode 100644 index 0000000000..550a951319 --- /dev/null +++ b/src/corecel/random/distribution/RoundedNonnegDistribution.hh @@ -0,0 +1,88 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file corecel/random/distribution/RoundedNonnegDistribution.hh +//! \sa RoundedNonnegDistribution.test.cc +//---------------------------------------------------------------------------// +#pragma once + +#include +#include + +#include "corecel/Macros.hh" +#include "corecel/Types.hh" +#include "corecel/math/Algorithms.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Sample from a distribution and return a rounded non-negative integer. + * + * This distribution wraps an arbitrary underlying sampler and rounds each + * value to the nearest integer by adding 0.5 and truncating. Negative values + * are clamped to zero, and values above the integer type range are clamped to + * the maximum representable value. + */ +template +class RoundedNonnegDistribution +{ + public: + //!@{ + //! \name Type aliases + using real_type = typename Distribution::real_type; + using result_type = IntType; + //!@} + + static_assert(std::is_floating_point_v); + static_assert(std::is_integral_v); + + public: + // Construct with distribution arguments + template + inline CELER_FUNCTION explicit RoundedNonnegDistribution(Args&&... args); + + // Sample a rounded non-negative integer + template + inline CELER_FUNCTION result_type operator()(Generator& rng); + + private: + Distribution sample_; +}; + +//---------------------------------------------------------------------------// +// INLINE DEFINITIONS +//---------------------------------------------------------------------------// +/*! + * Construct with arguments forwarded to the wrapped distribution. + */ +template +template +CELER_FUNCTION +RoundedNonnegDistribution::RoundedNonnegDistribution( + Args&&... args) + : sample_(celeritas::forward(args)...) +{ +} + +//---------------------------------------------------------------------------// +/*! + * Sample a random number according to the rounded non-negative distribution. + */ +template +template +CELER_FUNCTION auto +RoundedNonnegDistribution::operator()(Generator& rng) + -> result_type +{ + real_type value = sample_(rng) + real_type{0.5}; + value = clamp( + value, + real_type{0}, + static_cast(std::numeric_limits::max())); + return static_cast(value); +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/test/corecel/CMakeLists.txt b/test/corecel/CMakeLists.txt index 9f1c5d3dfd..78e83c8453 100644 --- a/test/corecel/CMakeLists.txt +++ b/test/corecel/CMakeLists.txt @@ -122,6 +122,7 @@ celeritas_add_test(random/distribution/NormalDistribution.test.cc) celeritas_add_test(random/distribution/PoissonDistribution.test.cc) celeritas_add_test(random/distribution/PowerDistribution.test.cc) celeritas_add_test(random/distribution/RadialDistribution.test.cc) +celeritas_add_test(random/distribution/RoundedNonnegDistribution.test.cc) celeritas_add_test(random/distribution/ReciprocalDistribution.test.cc) celeritas_add_test(random/distribution/RejectionSampler.test.cc) celeritas_add_test(random/distribution/TruncatedDistribution.test.cc) diff --git a/test/corecel/random/distribution/RoundedNonnegDistribution.test.cc b/test/corecel/random/distribution/RoundedNonnegDistribution.test.cc new file mode 100644 index 0000000000..bfdbc3b31f --- /dev/null +++ b/test/corecel/random/distribution/RoundedNonnegDistribution.test.cc @@ -0,0 +1,132 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file corecel/random/distribution/RoundedNonnegDistribution.test.cc +//---------------------------------------------------------------------------// +#include "corecel/random/distribution/RoundedNonnegDistribution.hh" + +#include +#include +#include +#include + +#include "corecel/cont/Range.hh" +#include "corecel/random/DiagnosticRngEngine.hh" +#include "corecel/random/Histogram.hh" +#include "corecel/random/distribution/UniformRealDistribution.hh" + +#include "celeritas_test.hh" + +namespace celeritas +{ +namespace test +{ +//---------------------------------------------------------------------------// + +namespace +{ +struct ConstantDistribution +{ + using real_type = double; + using result_type = double; + + double value; + + explicit CELER_FUNCTION ConstantDistribution(double v) : value(v) {} + + template + CELER_FUNCTION double operator()(Generator&) + { + return value; + } +}; + +struct CountingDistribution +{ + using real_type = double; + using result_type = double; + + int* count; + double value; + + CELER_FUNCTION CountingDistribution(int* c, double v) : count(c), value(v) + { + } + + template + CELER_FUNCTION double operator()(Generator&) + { + ++(*count); + return value; + } +}; +} // namespace + +TEST(RoundedNonnegDistributionTest, rounding) +{ + struct DummyRng + { + } rng; + + using RoundedConstant = RoundedNonnegDistribution; + static_assert( + std::is_same_v); + + EXPECT_EQ(0, RoundedConstant{-4.2}(rng)); + EXPECT_EQ(0, RoundedConstant{-0.2}(rng)); + EXPECT_EQ(0, RoundedConstant{0.49}(rng)); + EXPECT_EQ(1, RoundedConstant{0.5}(rng)); + EXPECT_EQ(2, RoundedConstant{1.5}(rng)); + EXPECT_EQ(3, RoundedConstant{2.5}(rng)); + + using RoundedU8 + = RoundedNonnegDistribution; + EXPECT_EQ(std::numeric_limits::max(), RoundedU8{1e6}(rng)); +} + +TEST(RoundedNonnegDistributionTest, one_sample_per_call) +{ + struct DummyRng + { + } rng; + + int num_samples = 10000; + int num_distribution_samples = 0; + RoundedNonnegDistribution sample_counting{ + &num_distribution_samples, 1.25}; + + for ([[maybe_unused]] int i : range(num_samples)) + { + EXPECT_EQ(1, sample_counting(rng)); + } + + EXPECT_EQ(num_samples, num_distribution_samples); +} + +TEST(RoundedNonnegDistributionTest, uniform) +{ + using RoundedUniform + = RoundedNonnegDistribution>; + + int num_samples = 10000; + DiagnosticRngEngine rng; + RoundedUniform sample_uniform{-0.8, 4.2}; + + Histogram histogram(5, {0, 5}); + for ([[maybe_unused]] int i : range(num_samples)) + { + histogram(static_cast(sample_uniform(rng))); + } + + static unsigned int const expected_counts[] + = {2673u, 1898u, 2078u, 1974u, 1377u}; + EXPECT_VEC_EQ(expected_counts, histogram.counts()); + EXPECT_EQ(0, histogram.underflow()); + EXPECT_EQ(0, histogram.overflow()); + EXPECT_EQ(2 * num_samples, rng.count()); +} + +//---------------------------------------------------------------------------// +} // namespace test +} // namespace celeritas From f9a689a12000c71c56bb5d59074e6f09b067af2d Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Fri, 22 May 2026 09:59:25 -0400 Subject: [PATCH 11/28] Update cherenkov baseline now that normal distribution can reuse values --- test/celeritas/optical/Cherenkov.test.cc | 34 +++++++++++------------- 1 file changed, 16 insertions(+), 18 deletions(-) diff --git a/test/celeritas/optical/Cherenkov.test.cc b/test/celeritas/optical/Cherenkov.test.cc index 8789f30a20..955ff3559b 100644 --- a/test/celeritas/optical/Cherenkov.test.cc +++ b/test/celeritas/optical/Cherenkov.test.cc @@ -439,13 +439,11 @@ TEST_F(CherenkovWaterTest, generator) // clang-format off static double const expected_costheta_dist[] - = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 52451, 10508, 0}; + = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 52173, 10417, 0}; static double const expected_energy_dist[] - = {3690, 3774, 3698, 3752, 3684, 3658, 3768, 3831, - 3921, 4029, 4025, 3941, 4134, 4286, 4307, 4461}; + = {3656, 3730, 3680, 3743, 3656, 3628, 3737, 3826, 3897, 3996, 4017, 3937, 4120, 4248, 4283, 4436}; static double const expected_displacement_dist[] - = {3909, 4064, 3802, 3920, 4001, 3904, 3891, 3955, - 3999, 3924, 3903, 3900, 3959, 3932, 4023, 3873}; + = {3869, 4048, 3777, 3886, 3970, 3898, 3862, 3944, 3965, 3893, 3899, 3863, 3937, 3911, 3997, 3871}; // clang-format on if (CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_DOUBLE) @@ -457,11 +455,11 @@ TEST_F(CherenkovWaterTest, generator) EXPECT_VEC_EQ(expected_costheta_dist, costheta_dist); EXPECT_VEC_EQ(expected_energy_dist, energy_dist); EXPECT_VEC_EQ(expected_displacement_dist, displacement_dist); - EXPECT_SOFT_EQ(0.73055857883146702, avg_costheta); - EXPECT_SOFT_EQ(4.0497726102182314e-06, avg_energy); - EXPECT_SOFT_EQ(0.50020101984474064, avg_displacement); - EXPECT_SOFT_EQ(983.734375, total_num_photons / num_samples); - EXPECT_SOFT_EQ(10.609603075017075, avg_engine_samples); + EXPECT_SOFT_EQ(0.73054872533349, avg_costheta); + EXPECT_SOFT_EQ(4.0511824907988e-06, avg_energy); + EXPECT_SOFT_EQ(0.50048096873882, avg_displacement); + EXPECT_SOFT_EQ(977.96875, total_num_photons / num_samples); + EXPECT_SOFT_EQ(10.607125738936, avg_engine_samples); } } @@ -483,11 +481,11 @@ TEST_F(CherenkovWaterTest, generator) Real3 pos = {sim.step_length(), 0, 0}; static double const expected_costheta_dist[] - = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 991}; + = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 997}; static double const expected_energy_dist[] - = {0, 0, 0, 0, 4, 14, 29, 26, 48, 51, 77, 103, 129, 132, 174, 204}; + = {0, 0, 0, 0, 5, 15, 30, 28, 52, 58, 78, 96, 125, 138, 165, 207}; static double const expected_displacement_dist[] = { - 123, 114, 103, 102, 83, 81, 80, 57, 60, 59, 31, 29, 36, 14, 16, 3}; + 121, 124, 94, 98, 80, 91, 80, 59, 59, 56, 33, 32, 37, 17, 14, 2}; sample(pre_step, particle, sim, pos, num_samples); @@ -496,11 +494,11 @@ TEST_F(CherenkovWaterTest, generator) EXPECT_VEC_EQ(expected_costheta_dist, costheta_dist); EXPECT_VEC_EQ(expected_energy_dist, energy_dist); EXPECT_VEC_EQ(expected_displacement_dist, displacement_dist); - EXPECT_SOFT_EQ(0.95045221539598979, avg_costheta); - EXPECT_SOFT_EQ(5.5902203966702514e-06, avg_energy); - EXPECT_SOFT_EQ(0.049715603846029896, avg_displacement); - EXPECT_SOFT_EQ(15.484375, total_num_photons / num_samples); - EXPECT_SOFT_EQ(25.077699293642784, avg_engine_samples); + EXPECT_SOFT_EQ(0.95082959158701, avg_costheta); + EXPECT_SOFT_EQ(5.5682878790974e-06, avg_energy); + EXPECT_SOFT_EQ(0.050055328486844, avg_displacement); + EXPECT_SOFT_EQ(15.578125, total_num_photons / num_samples); + EXPECT_SOFT_EQ(25.135406218656, avg_engine_samples); } } } From 83309538c29d285f5b441e2bdccd864273d4546f Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Fri, 22 May 2026 10:03:09 -0400 Subject: [PATCH 12/28] Use new distribution --- .../EnergyLossGaussianDistribution.hh | 1 - .../em/msc/detail/UrbanMscMinimalStepLimit.hh | 3 -- .../optical/gen/ScintillationGenerator.hh | 6 ++++ .../optical/gen/ScintillationOffload.hh | 7 +++-- .../distribution/PoissonDistribution.hh | 29 ++++++++++--------- .../distribution/RoundedNonnegDistribution.hh | 12 +++----- 6 files changed, 30 insertions(+), 28 deletions(-) diff --git a/src/celeritas/em/distribution/EnergyLossGaussianDistribution.hh b/src/celeritas/em/distribution/EnergyLossGaussianDistribution.hh index 9a3adab956..fab731a5c9 100644 --- a/src/celeritas/em/distribution/EnergyLossGaussianDistribution.hh +++ b/src/celeritas/em/distribution/EnergyLossGaussianDistribution.hh @@ -11,7 +11,6 @@ #include "corecel/Assert.hh" #include "corecel/Macros.hh" #include "corecel/Types.hh" -#include "corecel/math/Algorithms.hh" #include "corecel/random/distribution/NormalDistribution.hh" #include "corecel/random/distribution/TruncatedDistribution.hh" diff --git a/src/celeritas/em/msc/detail/UrbanMscMinimalStepLimit.hh b/src/celeritas/em/msc/detail/UrbanMscMinimalStepLimit.hh index fe2e45d752..3fad4f329e 100644 --- a/src/celeritas/em/msc/detail/UrbanMscMinimalStepLimit.hh +++ b/src/celeritas/em/msc/detail/UrbanMscMinimalStepLimit.hh @@ -6,13 +6,10 @@ //---------------------------------------------------------------------------// #pragma once -#include - #include "corecel/Macros.hh" #include "corecel/Types.hh" #include "corecel/math/Algorithms.hh" #include "corecel/random/distribution/NormalDistribution.hh" -#include "celeritas/Quantities.hh" #include "celeritas/Types.hh" #include "celeritas/em/data/UrbanMscData.hh" #include "celeritas/phys/PhysicsTrackView.hh" diff --git a/src/celeritas/optical/gen/ScintillationGenerator.hh b/src/celeritas/optical/gen/ScintillationGenerator.hh index 7f6ad91131..5734f14014 100644 --- a/src/celeritas/optical/gen/ScintillationGenerator.hh +++ b/src/celeritas/optical/gen/ScintillationGenerator.hh @@ -49,6 +49,12 @@ namespace optical * \note This performs the same sampling routine as in G4Scintillation class * of the Geant4 release 11.2 with some modifications. + * + * \note This could be made efficient by using separate generators + * for each scintillation component, and either using stratified sampling or + * poisson sampling to determine the number of photons for each component. It + * would also be a bit cleaner: separate generators would handle the energy + * distribution sampling. */ class ScintillationGenerator { diff --git a/src/celeritas/optical/gen/ScintillationOffload.hh b/src/celeritas/optical/gen/ScintillationOffload.hh index 46f39e3426..9a5abac523 100644 --- a/src/celeritas/optical/gen/ScintillationOffload.hh +++ b/src/celeritas/optical/gen/ScintillationOffload.hh @@ -10,6 +10,7 @@ #include "corecel/Macros.hh" #include "corecel/random/distribution/NormalDistribution.hh" #include "corecel/random/distribution/PoissonDistribution.hh" +#include "corecel/random/distribution/RoundedNonnegDistribution.hh" #include "celeritas/Quantities.hh" #include "celeritas/Types.hh" #include "celeritas/phys/ParticleTrackView.hh" @@ -131,9 +132,9 @@ ScintillationOffload::operator()(Generator& rng) real_type sigma = shared_.resolution_scale[pre_step_.material] * std::sqrt(mean_num_photons_); - result.num_photons = static_cast( - clamp_to_nonneg(NormalDistribution( - mean_num_photons_ + 0.5_r, sigma)(rng))); + result.num_photons + = RoundedNonnegDistribution>( + mean_num_photons_, sigma)(rng); } else if (mean_num_photons_ > 0) { diff --git a/src/corecel/random/distribution/PoissonDistribution.hh b/src/corecel/random/distribution/PoissonDistribution.hh index 168f381423..16659c642d 100644 --- a/src/corecel/random/distribution/PoissonDistribution.hh +++ b/src/corecel/random/distribution/PoissonDistribution.hh @@ -14,6 +14,7 @@ #include "corecel/math/Algorithms.hh" #include "NormalDistribution.hh" +#include "RoundedNonnegDistribution.hh" namespace celeritas { @@ -46,11 +47,11 @@ namespace celeritas * * In the degenerate case of \f$ \lambda = 0 \f$, the result is always zero. * - * \todo Rename to GeneralPoissonDistribution (or something similar) since + * \todo Rename to VariantPoissonDistribution (or something similar) since * it's effectively an inefficient variant combining: - * - an actual poisson distribution, - * - an "integer normal" distribution, and - * - a "zero" distribution. + * - an actual poisson distribution (using Knuth's method), + * - a rounded nonnegative distribution (for \f$ \lambda \gg 1 \f$), and + * - a delta distribution (for \f$ \lambda == 0 \f$ ). */ template class PoissonDistribution @@ -73,10 +74,13 @@ class PoissonDistribution template inline CELER_FUNCTION result_type operator()(Generator& rng); - //! Maximum value of lambda for using the direct method - static CELER_CONSTEXPR_FUNCTION int lambda_threshold() { return 16; } + //! Minimum value of lambda to approximate as a Gaussian + static CELER_CONSTEXPR_FUNCTION int large_lambda() { return 16; } private: + using RoundedNormal_t + = RoundedNonnegDistribution>; + enum class Method { zero, @@ -85,7 +89,7 @@ class PoissonDistribution }; Method method_; real_type exp_lambda_{}; - NormalDistribution sample_normal_{}; + RoundedNormal_t sample_normal_{}; }; //---------------------------------------------------------------------------// @@ -104,7 +108,7 @@ PoissonDistribution::PoissonDistribution(real_type lambda) { method_ = Method::zero; } - else if (lambda <= PoissonDistribution::lambda_threshold()) + else if (lambda <= PoissonDistribution::large_lambda()) { method_ = Method::poisson; exp_lambda_ = std::exp(lambda); @@ -112,8 +116,7 @@ PoissonDistribution::PoissonDistribution(real_type lambda) else { method_ = Method::gaussian; - // Add 0.5 to mean for correct rounding - sample_normal_ = NormalDistribution{lambda + 0.5, std::sqrt(lambda)}; + sample_normal_ = RoundedNormal_t{lambda, std::sqrt(lambda)}; } } @@ -141,10 +144,10 @@ CELER_FUNCTION auto PoissonDistribution::operator()(Generator& rng) return static_cast(k - 1); } case Method::gaussian: - // Use Gaussian approximation rounded to nearest integer - return static_cast( - clamp_to_nonneg(sample_normal_(rng))); + // Use Gaussian approximation rounded to nearest nonneg integer + return sample_normal_(rng); }; + CELER_ASSERT_UNREACHABLE(); } //---------------------------------------------------------------------------// diff --git a/src/corecel/random/distribution/RoundedNonnegDistribution.hh b/src/corecel/random/distribution/RoundedNonnegDistribution.hh index 550a951319..8a3846d073 100644 --- a/src/corecel/random/distribution/RoundedNonnegDistribution.hh +++ b/src/corecel/random/distribution/RoundedNonnegDistribution.hh @@ -3,7 +3,6 @@ // SPDX-License-Identifier: (Apache-2.0 OR MIT) //---------------------------------------------------------------------------// //! \file corecel/random/distribution/RoundedNonnegDistribution.hh -//! \sa RoundedNonnegDistribution.test.cc //---------------------------------------------------------------------------// #pragma once @@ -22,8 +21,9 @@ namespace celeritas * * This distribution wraps an arbitrary underlying sampler and rounds each * value to the nearest integer by adding 0.5 and truncating. Negative values - * are clamped to zero, and values above the integer type range are clamped to - * the maximum representable value. + * are clamped to zero. The underlying distribution MUST NOT return a value + * greater than the maximum representable value of the integer type (i.e., + * don't use this to sample the number of atoms in a gram of substance). */ template class RoundedNonnegDistribution @@ -76,11 +76,7 @@ CELER_FUNCTION auto RoundedNonnegDistribution::operator()(Generator& rng) -> result_type { - real_type value = sample_(rng) + real_type{0.5}; - value = clamp( - value, - real_type{0}, - static_cast(std::numeric_limits::max())); + real_type value = clamp_to_nonneg(sample_(rng) + real_type{0.5}); return static_cast(value); } From ef4a442c530c769c7043d5a09cd938de58b7646d Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Fri, 22 May 2026 10:47:19 -0400 Subject: [PATCH 13/28] Clearly differentiate non-scintillating materials in offload/inserter --- src/celeritas/optical/gen/ScintillationData.hh | 12 ++++++------ src/celeritas/optical/gen/ScintillationOffload.hh | 1 + .../optical/gen/detail/ScintSpectrumInserter.hh | 3 +++ 3 files changed, 10 insertions(+), 6 deletions(-) diff --git a/src/celeritas/optical/gen/ScintillationData.hh b/src/celeritas/optical/gen/ScintillationData.hh index 205269ebd2..23288b00d4 100644 --- a/src/celeritas/optical/gen/ScintillationData.hh +++ b/src/celeritas/optical/gen/ScintillationData.hh @@ -48,9 +48,10 @@ struct ScintDistributionRecord /*! * Unnormalized scintillation spectrum as a sum of independent components. * - * \todo The yield and resolution scale should live together (used for sampling - * the number of photons) and be separated from the normalized spectrum (used - * during generation, represented as a sum of components). + * \todo The yield and resolution scale should live together (used by \c + * ScintillationOffload for sampling the number of photons) and be separated + * from the normalized spectrum (used during generation by \c + * ScintillationGenerator, represented as a sum of components). * * - \c yield_per_energy is the average number of photons released by a unit of * locally deposited energy. @@ -64,11 +65,10 @@ struct ScintSpectrumRecord ItemRange yield_pdf; ItemRange components; - //! Whether all data are assigned and valid + //! Whether scintillation is present for this material explicit CELER_FUNCTION operator bool() const { - return yield_per_energy > 0 && !yield_pdf.empty() - && yield_pdf.size() == components.size(); + return yield_per_energy > 0; } }; diff --git a/src/celeritas/optical/gen/ScintillationOffload.hh b/src/celeritas/optical/gen/ScintillationOffload.hh index 9a5abac523..260babf110 100644 --- a/src/celeritas/optical/gen/ScintillationOffload.hh +++ b/src/celeritas/optical/gen/ScintillationOffload.hh @@ -109,6 +109,7 @@ CELER_FUNCTION ScintillationOffload::ScintillationOffload( //! \todo Use visible energy deposition when Birks law is implemented if (spectrum) { + // This material is a scintillator mean_num_photons_ = spectrum.yield_per_energy * energy_deposition.value(); } diff --git a/src/celeritas/optical/gen/detail/ScintSpectrumInserter.hh b/src/celeritas/optical/gen/detail/ScintSpectrumInserter.hh index 0a6c5d631e..c3ceae29bb 100644 --- a/src/celeritas/optical/gen/detail/ScintSpectrumInserter.hh +++ b/src/celeritas/optical/gen/detail/ScintSpectrumInserter.hh @@ -139,6 +139,9 @@ auto ScintSpectrumInserter::operator()(inp::ScintillationMaterial const& mat) spectrum.yield_per_energy = total_yield; spectrum.components = {begin_components, scint_records_.size_id()}; spectrum.yield_pdf = reals_.insert_back(yield_pdf.begin(), yield_pdf.end()); + CELER_ASSERT(total_yield > 0); + CELER_ASSERT(!spectrum.components.empty()); + CELER_ASSERT(spectrum.components.size() == spectrum.yield_pdf.size()); CELER_ENSURE(spectrum.components.size() == mat.components.size()); return spectra_.push_back(std::move(spectrum)); From 2703f70f1e01b8d28b2ca5261bce69098d1413f1 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Fri, 22 May 2026 10:51:55 -0400 Subject: [PATCH 14/28] Define separate knuth distribution --- .../EnergyLossUrbanDistribution.hh | 19 ++-- .../optical/gen/ScintillationOffload.hh | 20 ++-- .../random/distribution/NormalDistribution.hh | 15 +-- .../distribution/PoissonDistribution.hh | 107 ++++++++++++++---- .../distribution/PoissonDistribution.test.cc | 66 +++++++++++ 5 files changed, 180 insertions(+), 47 deletions(-) diff --git a/src/celeritas/em/distribution/EnergyLossUrbanDistribution.hh b/src/celeritas/em/distribution/EnergyLossUrbanDistribution.hh index 29527c2db9..e6f043ddd7 100644 --- a/src/celeritas/em/distribution/EnergyLossUrbanDistribution.hh +++ b/src/celeritas/em/distribution/EnergyLossUrbanDistribution.hh @@ -90,7 +90,10 @@ class EnergyLossUrbanDistribution static CELER_CONSTEXPR_FUNCTION real_type rate() { return 0.56; } //! Number of collisions above which to use faster sampling from Gaussian - static CELER_CONSTEXPR_FUNCTION size_type max_collisions() { return 8; } + static CELER_CONSTEXPR_FUNCTION size_type collision_poisson_threshold() + { + return 8; + } //! Threshold number of excitations used in width correction static CELER_CONSTEXPR_FUNCTION real_type exc_thresh() { return 42; } @@ -263,19 +266,20 @@ EnergyLossUrbanDistribution::sample_excitation_loss(Engine& rng) for (int i : range(2)) { - if (xs_exc_[i] > this->max_collisions()) + if (xs_exc_[i] > this->collision_poisson_threshold()) { // When the number of collisions is large, use faster approach // of sampling from a Gaussian mean += xs_exc_[i] * binding_energy_[i]; variance += xs_exc_[i] * ipow<2>(binding_energy_[i]); } - else + else if (xs_exc_[i] > 0) { // The loss due to excitation is \f$ \Delta E_{exc} = n_1 E_1 + n_2 // E_2 \f$, where the number of collisions \f$ n_i \f$ is sampled // from a Poisson distribution with mean \f$ \Sigma_i \f$ - unsigned int n = PoissonDistribution(xs_exc_[i])(rng); + unsigned int n + = PoissonDistributionKnuth(xs_exc_[i])(rng); if (n > 0) { UniformRealDistribution sample_fraction(n - 1, @@ -312,15 +316,16 @@ EnergyLossUrbanDistribution::sample_ionization_loss(Engine& rng) // Mean number of collisions in the fast simulation interval real_type mean_num_coll = 0; - if (xs_ion_ > this->max_collisions()) + if (xs_ion_ > this->collision_poisson_threshold()) { // When the number of collisions is large, fast sampling from a // Gaussian is used in the lower portion of the energy loss interval. // See PHYS332 section 2.4: Fast simulation for \f$ n_3 \ge 16 \f$ // Calculate the maximum value of \f$ \alpha \f$ (Eq. 25) - alpha = (xs_ion_ + this->max_collisions()) * energy_ratio - / (this->max_collisions() * energy_ratio + xs_ion_); + alpha + = (xs_ion_ + this->collision_poisson_threshold()) * energy_ratio + / (this->collision_poisson_threshold() * energy_ratio + xs_ion_); // Mean energy loss for a single collision of this type (Eq. 14) real_type const mean_loss_coll = alpha * std::log(alpha) / (alpha - 1); diff --git a/src/celeritas/optical/gen/ScintillationOffload.hh b/src/celeritas/optical/gen/ScintillationOffload.hh index 260babf110..d3340a062e 100644 --- a/src/celeritas/optical/gen/ScintillationOffload.hh +++ b/src/celeritas/optical/gen/ScintillationOffload.hh @@ -61,17 +61,15 @@ class ScintillationOffload private: units::ElementaryCharge charge_; - real_type step_length_; + real_type step_length_{}; OffloadPreStepData const& pre_step_; optical::GeneratorStepData post_step_; NativeCRef const& shared_; - real_type continuous_edep_fraction_; + real_type continuous_edep_fraction_{}; real_type mean_num_photons_{0}; - static CELER_CONSTEXPR_FUNCTION real_type poisson_threshold() - { - return 10; - } + // Use scaled gaussian above this average lambda + static constexpr real_type poisson_threshold = 10; }; //---------------------------------------------------------------------------// @@ -127,7 +125,7 @@ ScintillationOffload::operator()(Generator& rng) { // Material-only sampling optical::GeneratorDistributionData result; - if (mean_num_photons_ > poisson_threshold()) + if (mean_num_photons_ > poisson_threshold) { using namespace celeritas::literals; @@ -139,8 +137,12 @@ ScintillationOffload::operator()(Generator& rng) } else if (mean_num_photons_ > 0) { - result.num_photons = static_cast( - PoissonDistribution(mean_num_photons_)(rng)); + result.num_photons + = PoissonDistributionKnuth(mean_num_photons_)(rng); + } + else + { + result.num_photons = 0; } if (result.num_photons > 0) diff --git a/src/corecel/random/distribution/NormalDistribution.hh b/src/corecel/random/distribution/NormalDistribution.hh index 5a2ef371a3..93ad85436f 100644 --- a/src/corecel/random/distribution/NormalDistribution.hh +++ b/src/corecel/random/distribution/NormalDistribution.hh @@ -54,20 +54,17 @@ class NormalDistribution // Construct with mean and standard deviation inline CELER_FUNCTION NormalDistribution(real_type mean, real_type stddev); - //! Construct with unit deviation - explicit CELER_FUNCTION NormalDistribution(real_type mean) - : NormalDistribution{mean, 1} - { - } - // Construct from record explicit CELER_FUNCTION NormalDistribution(RecordT const& record) : NormalDistribution(record.mean, record.stddev) { } + //! Construct with unit deviation + explicit CELER_CEF NormalDistribution(real_type mean) : mean_{mean} {} + //! Construct with unit deviation and zero mean - CELER_FUNCTION NormalDistribution() : NormalDistribution{0, 1} {} + CELER_CEF NormalDistribution() = default; // Sample a random number according to the distribution template @@ -75,8 +72,8 @@ class NormalDistribution private: // Distribution properties - real_type mean_; - real_type stddev_; + real_type mean_{0}; + real_type stddev_{1}; // Intermediate samples real_type spare_{}; diff --git a/src/corecel/random/distribution/PoissonDistribution.hh b/src/corecel/random/distribution/PoissonDistribution.hh index 16659c642d..5069d30f1a 100644 --- a/src/corecel/random/distribution/PoissonDistribution.hh +++ b/src/corecel/random/distribution/PoissonDistribution.hh @@ -11,7 +11,6 @@ #include "corecel/Assert.hh" #include "corecel/Macros.hh" #include "corecel/Types.hh" -#include "corecel/math/Algorithms.hh" #include "NormalDistribution.hh" #include "RoundedNonnegDistribution.hh" @@ -20,7 +19,77 @@ namespace celeritas { //---------------------------------------------------------------------------// /*! - * Sample from a generalized Poisson distribution. + * Sample from a Poisson distribution using Knuth's algorithm. + * + * This algorithm should \em only be used for small, positive mean occurrences + * since the expected number of samples is proportional to the input value. + * + * See the \c PoissonDistribution below for documentation, as it should be used + * in the "general" case. + */ +template +class PoissonDistributionKnuth +{ + static_assert(std::is_floating_point_v); + + public: + //!@{ + //! \name Type aliases + using real_type = RealType; + using result_type = ::celeritas::size_type; + //!@} + public: + // Construct with distribution parameter + explicit inline CELER_FUNCTION PoissonDistributionKnuth(real_type lambda); + + //! Construct with default lambda of 1 + CELER_CEF PoissonDistributionKnuth() : exp_lambda_{constants::euler} {} + + // Sample a random number according to the distribution + template + inline CELER_FUNCTION result_type operator()(Generator& rng); + + private: + real_type exp_lambda_{}; +}; + +//---------------------------------------------------------------------------// +// INLINE DEFINITIONS +//---------------------------------------------------------------------------// +/*! + * Construct from the mean of the Poisson distribution. + */ +template +CELER_FUNCTION +PoissonDistributionKnuth::PoissonDistributionKnuth(real_type lambda) + : exp_lambda_{std::exp(lambda)} +{ + CELER_EXPECT(lambda >= 0); +} + +//---------------------------------------------------------------------------// +/*! + * Sample a random number according to the distribution. + */ +template +template +CELER_FUNCTION auto +PoissonDistributionKnuth::operator()(Generator& rng) -> result_type +{ + std::make_signed_t k{0}; + real_type p = exp_lambda_; + do + { + ++k; + p *= generate_canonical(rng); + } while (p > 1); + return static_cast(k - 1); + CELER_ASSERT_UNREACHABLE(); +} + +//---------------------------------------------------------------------------// +/*! + * Sample from a Poisson distribution. * * The Poisson distribution describes the probability of \f$ k \f$ events * occurring in a fixed interval given a mean rate of occurrence \f$ \lambda @@ -47,20 +116,21 @@ namespace celeritas * * In the degenerate case of \f$ \lambda = 0 \f$, the result is always zero. * - * \todo Rename to VariantPoissonDistribution (or something similar) since - * it's effectively an inefficient variant combining: + * \note This is effectivelty an inefficient variant combining: * - an actual poisson distribution (using Knuth's method), - * - a rounded nonnegative distribution (for \f$ \lambda \gg 1 \f$), and - * - a delta distribution (for \f$ \lambda == 0 \f$ ). + * - a rounded nonnegative normal distribution (for \f$ \lambda \gg 1 \f$), and + * - a delta distribution returning zero (for \f$ \lambda == 0 \f$ ). */ template class PoissonDistribution { + static_assert(std::is_floating_point_v); + public: //!@{ //! \name Type aliases using real_type = RealType; - using result_type = unsigned int; + using result_type = ::celeritas::size_type; //!@} public: @@ -78,17 +148,18 @@ class PoissonDistribution static CELER_CONSTEXPR_FUNCTION int large_lambda() { return 16; } private: + using PoissonKnuth_t = PoissonDistributionKnuth; using RoundedNormal_t - = RoundedNonnegDistribution>; + = RoundedNonnegDistribution, result_type>; enum class Method { zero, - poisson, + knuth, gaussian }; Method method_; - real_type exp_lambda_{}; + PoissonKnuth_t sample_knuth_{}; RoundedNormal_t sample_normal_{}; }; @@ -110,8 +181,8 @@ PoissonDistribution::PoissonDistribution(real_type lambda) } else if (lambda <= PoissonDistribution::large_lambda()) { - method_ = Method::poisson; - exp_lambda_ = std::exp(lambda); + method_ = Method::knuth; + sample_knuth_ = PoissonDistributionKnuth{lambda}; } else { @@ -133,16 +204,8 @@ CELER_FUNCTION auto PoissonDistribution::operator()(Generator& rng) { case Method::zero: return 0; - case Method::poisson: { - int k = 0; - real_type p = exp_lambda_; - do - { - ++k; - p *= generate_canonical(rng); - } while (p > 1); - return static_cast(k - 1); - } + case Method::knuth: + return sample_knuth_(rng); case Method::gaussian: // Use Gaussian approximation rounded to nearest nonneg integer return sample_normal_(rng); diff --git a/test/corecel/random/distribution/PoissonDistribution.test.cc b/test/corecel/random/distribution/PoissonDistribution.test.cc index b76b224b1b..7e80a3838a 100644 --- a/test/corecel/random/distribution/PoissonDistribution.test.cc +++ b/test/corecel/random/distribution/PoissonDistribution.test.cc @@ -9,7 +9,9 @@ #include "corecel/cont/Range.hh" #include "corecel/random/DiagnosticRngEngine.hh" #include "corecel/random/Histogram.hh" +#include "corecel/random/HistogramSampler.hh" +#include "TestMacros.hh" #include "celeritas_test.hh" namespace celeritas @@ -18,6 +20,70 @@ namespace test { //---------------------------------------------------------------------------// +TEST(PoissonDistributionKnuthTest, zero) +{ + HistogramSampler calc_histogram(8, {0, 8}, 100000); + + EXPECT_REF_EQ((SampledHistogram{{1, 0, 0, 0, 0, 0, 0, 0}, 2}), + calc_histogram(PoissonDistributionKnuth{0})); +} + +TEST(PoissonDistributionKnuthTest, small) +{ + HistogramSampler calc_histogram(8, {0, 8}, 100000); + std::vector actual; + + for (double lambda : {0.05, 0.1, 0.2, 0.5}) + { + PoissonDistributionKnuth sample_poisson{lambda}; + actual.push_back(calc_histogram(sample_poisson)); + } + + static SampledHistogram const expected[] = { + {{0.95199, 0.04675, 0.00124, 2e-05, 0, 0, 0, 0}, 2.09858}, + {{0.90498, 0.09051, 0.00437, 0.00014, 0, 0, 0, 0}, 2.19934}, + {{0.81752, 0.1651, 0.01628, 0.00104, 6e-05, 0, 0, 0}, 2.40204}, + {{0.60662, 0.30242, 0.07653, 0.01286, 0.0014, 0.00016, 1e-05, 0}, + 3.00104}, + }; + EXPECT_REF_EQ(expected, actual); +} + +TEST(PoissonDistributionKnuthTest, large) +{ + HistogramSampler calc_histogram(8, {0, 50}, 1000); + std::vector actual; + // Test default lambda=1 + actual.push_back(calc_histogram(PoissonDistributionKnuth{})); + + for (auto i : range(1, 18)) + { + PoissonDistributionKnuth sample_poisson{static_cast(i)}; + actual.push_back(calc_histogram(sample_poisson)); + } + static SampledHistogram const expected[] = { + {{0.15984, 0.00016, 0, 0, 0, 0, 0, 0}, 4.048}, + {{0.16, 0, 0, 0, 0, 0, 0, 0}, 4.036}, + {{0.15984, 0.00016, 0, 0, 0, 0, 0, 0}, 5.972}, + {{0.156, 0.004, 0, 0, 0, 0, 0, 0}, 7.74}, + {{0.14496, 0.01504, 0, 0, 0, 0, 0, 0}, 9.814}, + {{0.12064, 0.0392, 0.00016, 0, 0, 0, 0, 0}, 12.202}, + {{0.09984, 0.05856, 0.0016, 0, 0, 0, 0, 0}, 13.876}, + {{0.07392, 0.08288, 0.0032, 0, 0, 0, 0, 0}, 15.78}, + {{0.04768, 0.10112, 0.01088, 0.00032, 0, 0, 0, 0}, 18.258}, + {{0.03344, 0.10608, 0.02032, 0.00016, 0, 0, 0, 0}, 20.07}, + {{0.02144, 0.10528, 0.03248, 0.0008, 0, 0, 0, 0}, 21.976}, + {{0.01264, 0.1024, 0.04272, 0.00208, 0.00016, 0, 0, 0}, 23.454}, + {{0.00944, 0.08576, 0.05984, 0.00496, 0, 0, 0, 0}, 25.636}, + {{0.00336, 0.0696, 0.0752, 0.0112, 0.00064, 0, 0, 0}, 28.374}, + {{0.00208, 0.05632, 0.08352, 0.01744, 0.00064, 0, 0, 0}, 29.918}, + {{0.00128, 0.04288, 0.0824, 0.03104, 0.00224, 0.00016, 0, 0}, 32.384}, + {{0.0008, 0.0288, 0.08672, 0.03984, 0.00384, 0, 0, 0}, 34.208}, + {{0.00032, 0.02192, 0.08336, 0.04736, 0.00672, 0.00032, 0, 0}, 36.12}, + }; + EXPECT_REF_EQ(expected, actual) << repr(actual); +} + TEST(PoissonDistributionTest, bin_zero) { int num_samples = 100; From 7bb01bab65e910616d553d447a327d8a361aa364 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Fri, 22 May 2026 11:29:15 -0400 Subject: [PATCH 15/28] Address feedback --- .../em/distribution/EnergyLossUrbanDistribution.hh | 3 +-- .../random/distribution/PoissonDistribution.hh | 11 ++++++----- .../random/distribution/TruncatedDistribution.hh | 12 ++++++------ .../distribution/TruncatedDistribution.test.cc | 2 +- 4 files changed, 14 insertions(+), 14 deletions(-) diff --git a/src/celeritas/em/distribution/EnergyLossUrbanDistribution.hh b/src/celeritas/em/distribution/EnergyLossUrbanDistribution.hh index e6f043ddd7..dd618304d0 100644 --- a/src/celeritas/em/distribution/EnergyLossUrbanDistribution.hh +++ b/src/celeritas/em/distribution/EnergyLossUrbanDistribution.hh @@ -278,8 +278,7 @@ EnergyLossUrbanDistribution::sample_excitation_loss(Engine& rng) // The loss due to excitation is \f$ \Delta E_{exc} = n_1 E_1 + n_2 // E_2 \f$, where the number of collisions \f$ n_i \f$ is sampled // from a Poisson distribution with mean \f$ \Sigma_i \f$ - unsigned int n - = PoissonDistributionKnuth(xs_exc_[i])(rng); + auto n = PoissonDistributionKnuth(xs_exc_[i])(rng); if (n > 0) { UniformRealDistribution sample_fraction(n - 1, diff --git a/src/corecel/random/distribution/PoissonDistribution.hh b/src/corecel/random/distribution/PoissonDistribution.hh index 5069d30f1a..ef8a786d39 100644 --- a/src/corecel/random/distribution/PoissonDistribution.hh +++ b/src/corecel/random/distribution/PoissonDistribution.hh @@ -9,9 +9,11 @@ #include #include "corecel/Assert.hh" +#include "corecel/Constants.hh" #include "corecel/Macros.hh" #include "corecel/Types.hh" +#include "GenerateCanonical.hh" #include "NormalDistribution.hh" #include "RoundedNonnegDistribution.hh" @@ -84,7 +86,6 @@ PoissonDistributionKnuth::operator()(Generator& rng) -> result_type p *= generate_canonical(rng); } while (p > 1); return static_cast(k - 1); - CELER_ASSERT_UNREACHABLE(); } //---------------------------------------------------------------------------// @@ -116,7 +117,7 @@ PoissonDistributionKnuth::operator()(Generator& rng) -> result_type * * In the degenerate case of \f$ \lambda = 0 \f$, the result is always zero. * - * \note This is effectivelty an inefficient variant combining: + * \note This is effectively a rough-and-ready variant selecting between: * - an actual poisson distribution (using Knuth's method), * - a rounded nonnegative normal distribution (for \f$ \lambda \gg 1 \f$), and * - a delta distribution returning zero (for \f$ \lambda == 0 \f$ ). @@ -156,7 +157,7 @@ class PoissonDistribution { zero, knuth, - gaussian + normal }; Method method_; PoissonKnuth_t sample_knuth_{}; @@ -186,7 +187,7 @@ PoissonDistribution::PoissonDistribution(real_type lambda) } else { - method_ = Method::gaussian; + method_ = Method::normal; sample_normal_ = RoundedNormal_t{lambda, std::sqrt(lambda)}; } } @@ -206,7 +207,7 @@ CELER_FUNCTION auto PoissonDistribution::operator()(Generator& rng) return 0; case Method::knuth: return sample_knuth_(rng); - case Method::gaussian: + case Method::normal: // Use Gaussian approximation rounded to nearest nonneg integer return sample_normal_(rng); }; diff --git a/src/corecel/random/distribution/TruncatedDistribution.hh b/src/corecel/random/distribution/TruncatedDistribution.hh index 40a8cb7fa3..19c61c3343 100644 --- a/src/corecel/random/distribution/TruncatedDistribution.hh +++ b/src/corecel/random/distribution/TruncatedDistribution.hh @@ -28,7 +28,7 @@ namespace celeritas * prevents more than 32 failed samples. Since GPU performance is so sensitive * to rejection failures, consider transforming the underlying distribution if * this limit is hit. When using this distribution for physics distributions, - * it is is \em strongly recommended to use \c test::HistogramSampler to verify + * it is \em strongly recommended to use \c test::HistogramSampler to verify * the number of samples being taken. */ template @@ -94,22 +94,22 @@ template CELER_FUNCTION auto TruncatedDistribution::operator()(Generator& rng) -> result_type { - int num_remaining_samples = max_debug_samples; + int num_remaining_samples = max_debug_samples + 1; result_type result; do { - // Reject samples outside the truncation bounds - result = sample_(rng); - // Prevent infinite loops (debug assertions only) if constexpr (CELERITAS_DEBUG) { - if (--num_remaining_samples < 0) + // Prevent infinite loops (debug assertions only) + if (--num_remaining_samples == 0) { CELER_DEBUG_FAIL( "too many samples taken in TruncatedDistribution", internal); } } + + result = sample_(rng); } while (result < lower_ || result > upper_); return result; } diff --git a/test/corecel/random/distribution/TruncatedDistribution.test.cc b/test/corecel/random/distribution/TruncatedDistribution.test.cc index 80231c0bfd..0008307c2e 100644 --- a/test/corecel/random/distribution/TruncatedDistribution.test.cc +++ b/test/corecel/random/distribution/TruncatedDistribution.test.cc @@ -75,7 +75,7 @@ TEST(TruncatedDistributionTest, TEST_IF_CELERITAS_DEBUG(rejection)) DiagnosticRngEngine rng; EXPECT_THROW(sample_forever(rng), DebugError); - EXPECT_EQ(66, rng.exchange_count()); + EXPECT_EQ(64, rng.exchange_count()); } //---------------------------------------------------------------------------// From 0a31daaf1722cfa65a7421ccc3d8a72cf44a3a45 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Fri, 22 May 2026 13:52:35 -0400 Subject: [PATCH 16/28] Eliminate unnecessary real type and use cold code for assertion checking --- src/corecel/random/distribution/DeltaDistribution.hh | 3 +-- .../random/distribution/DistributionTypeTraits.hh | 1 - src/corecel/random/distribution/TruncatedDistribution.hh | 9 ++++----- 3 files changed, 5 insertions(+), 8 deletions(-) diff --git a/src/corecel/random/distribution/DeltaDistribution.hh b/src/corecel/random/distribution/DeltaDistribution.hh index ad76236da7..54c2a9256f 100644 --- a/src/corecel/random/distribution/DeltaDistribution.hh +++ b/src/corecel/random/distribution/DeltaDistribution.hh @@ -6,7 +6,6 @@ //---------------------------------------------------------------------------// #pragma once -#include "corecel/Types.hh" #include "corecel/random/data/DistributionData.hh" namespace celeritas @@ -31,7 +30,7 @@ class DeltaDistribution { } - // Construct with record + //! Construct from a record explicit CELER_FUNCTION DeltaDistribution(RecordT const& record) : value_(record.value) { diff --git a/src/corecel/random/distribution/DistributionTypeTraits.hh b/src/corecel/random/distribution/DistributionTypeTraits.hh index cbc82ad957..f8b8205ada 100644 --- a/src/corecel/random/distribution/DistributionTypeTraits.hh +++ b/src/corecel/random/distribution/DistributionTypeTraits.hh @@ -9,7 +9,6 @@ #include "corecel/Types.hh" #include "corecel/cont/EnumClassUtils.hh" #include "corecel/random/Types.hh" -#include "corecel/random/data/DistributionData.hh" #include "DeltaDistribution.hh" #include "IsotropicDistribution.hh" diff --git a/src/corecel/random/distribution/TruncatedDistribution.hh b/src/corecel/random/distribution/TruncatedDistribution.hh index 19c61c3343..9fa85fad1b 100644 --- a/src/corecel/random/distribution/TruncatedDistribution.hh +++ b/src/corecel/random/distribution/TruncatedDistribution.hh @@ -37,7 +37,6 @@ class TruncatedDistribution public: //!@{ //! \name Type aliases - using real_type = typename Distribution::real_type; using result_type = typename Distribution::result_type; using RecordT = TruncatedDistributionRecord; //!@} @@ -48,7 +47,7 @@ class TruncatedDistribution inline CELER_FUNCTION TruncatedDistribution(T lower, U upper, Args&&... args); - // Construct from record + //! Construct from a device-friendly variant record explicit CELER_FUNCTION TruncatedDistribution(RecordT const& record) : TruncatedDistribution( record.lower, record.upper, Distribution(record.distribution)) @@ -64,8 +63,8 @@ class TruncatedDistribution private: Distribution sample_; - real_type lower_; - real_type upper_; + result_type lower_; + result_type upper_; }; //---------------------------------------------------------------------------// @@ -101,7 +100,7 @@ TruncatedDistribution::operator()(Generator& rng) -> result_type if constexpr (CELERITAS_DEBUG) { // Prevent infinite loops (debug assertions only) - if (--num_remaining_samples == 0) + if (CELER_UNLIKELY(--num_remaining_samples == 0)) { CELER_DEBUG_FAIL( "too many samples taken in TruncatedDistribution", From 4e71c80e71988548981be0d25b3fbfc8ee113a29 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Fri, 22 May 2026 14:42:48 -0400 Subject: [PATCH 17/28] Address feedback --- .../optical/gen/ScintillationGenerator.hh | 6 ------ .../distribution/PoissonDistribution.hh | 20 ++++++++++++++----- 2 files changed, 15 insertions(+), 11 deletions(-) diff --git a/src/celeritas/optical/gen/ScintillationGenerator.hh b/src/celeritas/optical/gen/ScintillationGenerator.hh index 5734f14014..7f6ad91131 100644 --- a/src/celeritas/optical/gen/ScintillationGenerator.hh +++ b/src/celeritas/optical/gen/ScintillationGenerator.hh @@ -49,12 +49,6 @@ namespace optical * \note This performs the same sampling routine as in G4Scintillation class * of the Geant4 release 11.2 with some modifications. - * - * \note This could be made efficient by using separate generators - * for each scintillation component, and either using stratified sampling or - * poisson sampling to determine the number of photons for each component. It - * would also be a bit cleaner: separate generators would handle the energy - * distribution sampling. */ class ScintillationGenerator { diff --git a/src/corecel/random/distribution/PoissonDistribution.hh b/src/corecel/random/distribution/PoissonDistribution.hh index ef8a786d39..0a648f9c83 100644 --- a/src/corecel/random/distribution/PoissonDistribution.hh +++ b/src/corecel/random/distribution/PoissonDistribution.hh @@ -21,13 +21,18 @@ namespace celeritas { //---------------------------------------------------------------------------// /*! - * Sample from a Poisson distribution using Knuth's algorithm. + * Sample from a Poisson distribution *for small N* using Knuth's algorithm. * * This algorithm should \em only be used for small, positive mean occurrences * since the expected number of samples is proportional to the input value. * - * See the \c PoissonDistribution below for documentation, as it should be used - * in the "general" case. + * A hard-coded maximum limit prevents developers from unwittingly use this for + * large values. If you hit the related assertion, use \c PoissonDistribution + * below or a different algorithm with a more accurate approximation (e.g., + * Wilson-Hilferty). + * + * See the \c PoissonDistribution below for documentation of the Poisson + * distribution, as it should be used in the "general" case. */ template class PoissonDistributionKnuth @@ -51,6 +56,9 @@ class PoissonDistributionKnuth template inline CELER_FUNCTION result_type operator()(Generator& rng); + //! Maximum limit to prevent developers from getting into trouble + static constexpr real_type too_expensive_lambda{32}; + private: real_type exp_lambda_{}; }; @@ -67,6 +75,7 @@ PoissonDistributionKnuth::PoissonDistributionKnuth(real_type lambda) : exp_lambda_{std::exp(lambda)} { CELER_EXPECT(lambda >= 0); + CELER_EXPECT(lambda < too_expensive_lambda); // See class docs } //---------------------------------------------------------------------------// @@ -112,10 +121,11 @@ PoissonDistributionKnuth::operator()(Generator& rng) -> result_type * used for large \f$ \lambda \f$. * * Geant4 uses Knuth's algorithm for \f$ \lambda \le 16 \f$ and a Gaussian - * approximation for \f$ \lambda > 16 \f$ (see \c G4Poisson), which is faster + * approximation for \f$ \lambda > 16 \f$ (see \c G4Poisson ), which is faster * but less accurate than other methods. The same approach is used here. * - * In the degenerate case of \f$ \lambda = 0 \f$, the result is always zero. + * In the degenerate case of \f$ \lambda = 0 \f$, the result is always zero and + * requires no random numbers to be drawn. * * \note This is effectively a rough-and-ready variant selecting between: * - an actual poisson distribution (using Knuth's method), From 4c0c6ed493f23fbf095f068e91dab34c96a7ad48 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Sat, 23 May 2026 07:43:21 -0400 Subject: [PATCH 18/28] Fix implicit casting and add max assertion --- .../random/distribution/PoissonDistribution.hh | 1 + .../distribution/RoundedNonnegDistribution.hh | 3 ++- .../distribution/PoissonDistribution.test.cc | 16 ++++++++++++---- 3 files changed, 15 insertions(+), 5 deletions(-) diff --git a/src/corecel/random/distribution/PoissonDistribution.hh b/src/corecel/random/distribution/PoissonDistribution.hh index 0a648f9c83..a35dfe68e2 100644 --- a/src/corecel/random/distribution/PoissonDistribution.hh +++ b/src/corecel/random/distribution/PoissonDistribution.hh @@ -45,6 +45,7 @@ class PoissonDistributionKnuth using real_type = RealType; using result_type = ::celeritas::size_type; //!@} + public: // Construct with distribution parameter explicit inline CELER_FUNCTION PoissonDistributionKnuth(real_type lambda); diff --git a/src/corecel/random/distribution/RoundedNonnegDistribution.hh b/src/corecel/random/distribution/RoundedNonnegDistribution.hh index 8a3846d073..57f4830e78 100644 --- a/src/corecel/random/distribution/RoundedNonnegDistribution.hh +++ b/src/corecel/random/distribution/RoundedNonnegDistribution.hh @@ -6,12 +6,12 @@ //---------------------------------------------------------------------------// #pragma once -#include #include #include "corecel/Macros.hh" #include "corecel/Types.hh" #include "corecel/math/Algorithms.hh" +#include "corecel/math/NumericLimits.hh" namespace celeritas { @@ -77,6 +77,7 @@ RoundedNonnegDistribution::operator()(Generator& rng) -> result_type { real_type value = clamp_to_nonneg(sample_(rng) + real_type{0.5}); + CELER_ENSURE(value < static_cast(NumericLimits::max())); return static_cast(value); } diff --git a/test/corecel/random/distribution/PoissonDistribution.test.cc b/test/corecel/random/distribution/PoissonDistribution.test.cc index 7e80a3838a..2270e84c6d 100644 --- a/test/corecel/random/distribution/PoissonDistribution.test.cc +++ b/test/corecel/random/distribution/PoissonDistribution.test.cc @@ -18,6 +18,13 @@ namespace celeritas { namespace test { +//---------------------------------------------------------------------------// +// SampleHistogram bins doubles, but Knuth poisson returns integers +CELER_FORCEINLINE double static_cast_double(size_type v) +{ + return static_cast(v); +}; + //---------------------------------------------------------------------------// TEST(PoissonDistributionKnuthTest, zero) @@ -25,7 +32,8 @@ TEST(PoissonDistributionKnuthTest, zero) HistogramSampler calc_histogram(8, {0, 8}, 100000); EXPECT_REF_EQ((SampledHistogram{{1, 0, 0, 0, 0, 0, 0, 0}, 2}), - calc_histogram(PoissonDistributionKnuth{0})); + calc_histogram(static_cast_double, + PoissonDistributionKnuth{0})); } TEST(PoissonDistributionKnuthTest, small) @@ -36,7 +44,7 @@ TEST(PoissonDistributionKnuthTest, small) for (double lambda : {0.05, 0.1, 0.2, 0.5}) { PoissonDistributionKnuth sample_poisson{lambda}; - actual.push_back(calc_histogram(sample_poisson)); + actual.push_back(calc_histogram(static_cast_double, sample_poisson)); } static SampledHistogram const expected[] = { @@ -51,7 +59,7 @@ TEST(PoissonDistributionKnuthTest, small) TEST(PoissonDistributionKnuthTest, large) { - HistogramSampler calc_histogram(8, {0, 50}, 1000); + HistogramSampler calc_histogram(8, {0, 50}, 10000); std::vector actual; // Test default lambda=1 actual.push_back(calc_histogram(PoissonDistributionKnuth{})); @@ -59,7 +67,7 @@ TEST(PoissonDistributionKnuthTest, large) for (auto i : range(1, 18)) { PoissonDistributionKnuth sample_poisson{static_cast(i)}; - actual.push_back(calc_histogram(sample_poisson)); + actual.push_back(calc_histogram(static_cast_double, sample_poisson)); } static SampledHistogram const expected[] = { {{0.15984, 0.00016, 0, 0, 0, 0, 0, 0}, 4.048}, From e0dcfa4cc7240259b43d58f0f8e2bb3899e62fcc Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Thu, 28 May 2026 07:46:09 -0400 Subject: [PATCH 19/28] Add missing numeric limits --- src/corecel/math/NumericLimits.hh | 32 +++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/src/corecel/math/NumericLimits.hh b/src/corecel/math/NumericLimits.hh index c8d5900382..77ee8e00b4 100644 --- a/src/corecel/math/NumericLimits.hh +++ b/src/corecel/math/NumericLimits.hh @@ -58,6 +58,22 @@ struct numeric_limits SCCEF_ double infinity() { return __builtin_huge_val(); } }; +template<> +struct numeric_limits +{ + SCCEF_ signed char lowest() { return CHAR_MIN; } + SCCEF_ signed char min() { return CHAR_MIN; } + SCCEF_ signed char max() { return CHAR_MAX; } +}; + +template<> +struct numeric_limits +{ + SCCEF_ short lowest() { return SHRT_MIN; } + SCCEF_ short min() { return SHRT_MIN; } + SCCEF_ short max() { return SHRT_MAX; } +}; + template<> struct numeric_limits { @@ -82,6 +98,22 @@ struct numeric_limits SCCEF_ long long max() { return LLONG_MAX; } }; +template<> +struct numeric_limits +{ + SCCEF_ unsigned char lowest() { return 0; } + SCCEF_ unsigned char min() { return 0; } + SCCEF_ unsigned char max() { return UCHAR_MAX; } +}; + +template<> +struct numeric_limits +{ + SCCEF_ unsigned short lowest() { return 0; } + SCCEF_ unsigned short min() { return 0; } + SCCEF_ unsigned short max() { return USHRT_MAX; } +}; + template<> struct numeric_limits { From e7a41df0d60a2abb417703ca4f57fa32f3f34ea7 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Thu, 28 May 2026 11:15:48 -0400 Subject: [PATCH 20/28] Fix int types --- src/corecel/math/NumericLimits.hh | 22 ++++------------------ 1 file changed, 4 insertions(+), 18 deletions(-) diff --git a/src/corecel/math/NumericLimits.hh b/src/corecel/math/NumericLimits.hh index 77ee8e00b4..a97d80ea36 100644 --- a/src/corecel/math/NumericLimits.hh +++ b/src/corecel/math/NumericLimits.hh @@ -26,6 +26,8 @@ namespace celeritas * \deprecated The name of this class will change to NumericLimits to conform * to the Celeritas naming system, since it is not a complete replacement for * the \c std implementation. + * + * \todo Replace with the cuda/hip standard libraries or a modified llvm port */ template struct numeric_limits; @@ -58,22 +60,6 @@ struct numeric_limits SCCEF_ double infinity() { return __builtin_huge_val(); } }; -template<> -struct numeric_limits -{ - SCCEF_ signed char lowest() { return CHAR_MIN; } - SCCEF_ signed char min() { return CHAR_MIN; } - SCCEF_ signed char max() { return CHAR_MAX; } -}; - -template<> -struct numeric_limits -{ - SCCEF_ short lowest() { return SHRT_MIN; } - SCCEF_ short min() { return SHRT_MIN; } - SCCEF_ short max() { return SHRT_MAX; } -}; - template<> struct numeric_limits { @@ -103,7 +89,7 @@ struct numeric_limits { SCCEF_ unsigned char lowest() { return 0; } SCCEF_ unsigned char min() { return 0; } - SCCEF_ unsigned char max() { return UCHAR_MAX; } + SCCEF_ unsigned char max() { return (unsigned char)(UCHAR_MAX); } }; template<> @@ -111,7 +97,7 @@ struct numeric_limits { SCCEF_ unsigned short lowest() { return 0; } SCCEF_ unsigned short min() { return 0; } - SCCEF_ unsigned short max() { return USHRT_MAX; } + SCCEF_ unsigned short max() { return (unsigned short)(USHRT_MAX); } }; template<> From 8c2d5cabd017114ffd95798f63d309ef7e80996f Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Thu, 28 May 2026 18:26:13 -0400 Subject: [PATCH 21/28] Fix tests --- .../random/distribution/PoissonDistribution.test.cc | 2 +- .../distribution/RoundedNonnegDistribution.test.cc | 9 ++++++--- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/test/corecel/random/distribution/PoissonDistribution.test.cc b/test/corecel/random/distribution/PoissonDistribution.test.cc index 2270e84c6d..12f8381cea 100644 --- a/test/corecel/random/distribution/PoissonDistribution.test.cc +++ b/test/corecel/random/distribution/PoissonDistribution.test.cc @@ -59,7 +59,7 @@ TEST(PoissonDistributionKnuthTest, small) TEST(PoissonDistributionKnuthTest, large) { - HistogramSampler calc_histogram(8, {0, 50}, 10000); + HistogramSampler calc_histogram(8, {0, 50}, 1000); std::vector actual; // Test default lambda=1 actual.push_back(calc_histogram(PoissonDistributionKnuth{})); diff --git a/test/corecel/random/distribution/RoundedNonnegDistribution.test.cc b/test/corecel/random/distribution/RoundedNonnegDistribution.test.cc index bfdbc3b31f..02804e9dd6 100644 --- a/test/corecel/random/distribution/RoundedNonnegDistribution.test.cc +++ b/test/corecel/random/distribution/RoundedNonnegDistribution.test.cc @@ -80,9 +80,12 @@ TEST(RoundedNonnegDistributionTest, rounding) EXPECT_EQ(2, RoundedConstant{1.5}(rng)); EXPECT_EQ(3, RoundedConstant{2.5}(rng)); - using RoundedU8 - = RoundedNonnegDistribution; - EXPECT_EQ(std::numeric_limits::max(), RoundedU8{1e6}(rng)); + if (CELERITAS_DEBUG) + { + using RoundedU8 + = RoundedNonnegDistribution; + EXPECT_THROW(RoundedU8{1e6}(rng), DebugError); + } } TEST(RoundedNonnegDistributionTest, one_sample_per_call) From c434dc4f6f552222df32358f9dd32c31bbd9174a Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Fri, 29 May 2026 14:32:26 -0400 Subject: [PATCH 22/28] Fix windows error --- src/celeritas/em/distribution/EnergyLossUrbanDistribution.hh | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/celeritas/em/distribution/EnergyLossUrbanDistribution.hh b/src/celeritas/em/distribution/EnergyLossUrbanDistribution.hh index dd618304d0..941082289b 100644 --- a/src/celeritas/em/distribution/EnergyLossUrbanDistribution.hh +++ b/src/celeritas/em/distribution/EnergyLossUrbanDistribution.hh @@ -281,9 +281,8 @@ EnergyLossUrbanDistribution::sample_excitation_loss(Engine& rng) auto n = PoissonDistributionKnuth(xs_exc_[i])(rng); if (n > 0) { - UniformRealDistribution sample_fraction(n - 1, - n + 1); - result += sample_fraction(rng) * binding_energy_[i]; + UniformRealDistribution sample_smooth(-1, 1); + result += (n + sample_smooth(rng)) * binding_energy_[i]; } } } From 3015e2fc2335a1c6a057dac219fb82a554333965 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Wed, 3 Jun 2026 11:19:14 -0400 Subject: [PATCH 23/28] one more static cast --- .../random/distribution/PoissonDistribution.test.cc | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/corecel/random/distribution/PoissonDistribution.test.cc b/test/corecel/random/distribution/PoissonDistribution.test.cc index 12f8381cea..3764e85bff 100644 --- a/test/corecel/random/distribution/PoissonDistribution.test.cc +++ b/test/corecel/random/distribution/PoissonDistribution.test.cc @@ -23,7 +23,7 @@ namespace test CELER_FORCEINLINE double static_cast_double(size_type v) { return static_cast(v); -}; +} //---------------------------------------------------------------------------// @@ -104,7 +104,7 @@ TEST(PoissonDistributionTest, bin_zero) Histogram histogram(4, {0, 1e-3}); for ([[maybe_unused]] int i : range(num_samples)) { - histogram(sample_poisson(rng)); + histogram(static_cast_double(sample_poisson(rng))); } static unsigned int const expected_counts[] = {100, 0, 0, 0}; EXPECT_VEC_EQ(expected_counts, histogram.counts()); @@ -124,7 +124,7 @@ TEST(PoissonDistributionTest, bin_small) Histogram histogram(16, {0, 16}); for ([[maybe_unused]] int i : range(num_samples)) { - histogram(sample_poisson(rng)); + histogram(static_cast_double(sample_poisson(rng))); } static unsigned int const expected_counts[] = { 177, 762, 1444, 1971, 1950, 1586, 1054, 562, 286, 125, 55, 18, 5, 1, 3, 1}; @@ -145,7 +145,7 @@ TEST(PoissonDistributionTest, bin_large) Histogram histogram(60, {34.5, 94.5}); for ([[maybe_unused]] int i : range(num_samples)) { - histogram(static_cast(sample_poisson(rng))); + histogram(static_cast_double(sample_poisson(rng))); } static unsigned int const expected_counts[] = {1, 1, 5, 2, 5, 6, 6, 11, 11, 11, 28, 45, From 486a2269811975ebe5b2b5f7a3bc300e6b4f61d5 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Thu, 4 Jun 2026 09:17:58 -0400 Subject: [PATCH 24/28] Explicitly cast from int to double --- .../em/distribution/EnergyLossUrbanDistribution.hh | 3 ++- src/celeritas/ext/detail/GeantProcessImporter.cc | 11 ++++++----- src/corecel/cont/Span.hh | 5 ++--- src/corecel/grid/UniformGrid.hh | 2 +- src/corecel/grid/UniformGridData.hh | 3 ++- src/corecel/grid/VectorUtils.hh | 3 ++- src/geocel/rasterize/SafetyImager.t.hh | 5 +++-- test/corecel/random/Histogram.hh | 3 ++- .../random/distribution/PoissonDistribution.test.cc | 3 ++- test/testdetail/TestMacrosImpl.cc | 5 ++--- test/testdetail/TestMacrosImpl.hh | 7 +++---- 11 files changed, 27 insertions(+), 23 deletions(-) diff --git a/src/celeritas/em/distribution/EnergyLossUrbanDistribution.hh b/src/celeritas/em/distribution/EnergyLossUrbanDistribution.hh index 941082289b..7ed4c922e9 100644 --- a/src/celeritas/em/distribution/EnergyLossUrbanDistribution.hh +++ b/src/celeritas/em/distribution/EnergyLossUrbanDistribution.hh @@ -282,7 +282,8 @@ EnergyLossUrbanDistribution::sample_excitation_loss(Engine& rng) if (n > 0) { UniformRealDistribution sample_smooth(-1, 1); - result += (n + sample_smooth(rng)) * binding_energy_[i]; + result += (static_cast(n) + sample_smooth(rng)) + * binding_energy_[i]; } } } diff --git a/src/celeritas/ext/detail/GeantProcessImporter.cc b/src/celeritas/ext/detail/GeantProcessImporter.cc index d1e613adce..674c6a22f9 100644 --- a/src/celeritas/ext/detail/GeantProcessImporter.cc +++ b/src/celeritas/ext/detail/GeantProcessImporter.cc @@ -35,7 +35,9 @@ #include "corecel/Assert.hh" #include "corecel/cont/Range.hh" +#include "corecel/cont/Span.hh" #include "corecel/data/HyperslabIndexer.hh" +#include "corecel/grid/VectorUtils.hh" #include "corecel/io/Logger.hh" #include "corecel/math/Algorithms.hh" #include "corecel/math/SoftEqual.hh" @@ -387,19 +389,18 @@ inp::UniformGrid import_physics_log_vector(G4PhysicsVector const& pv, double const x_scaling = native_value_from_clhep(units[0]); double const y_scaling = native_value_from_clhep(units[1]); auto size = pv.GetVectorLength(); + double const bounds[] = {pv.Energy(0), pv.Energy(size - 1)}; inp::UniformGrid grid; - grid.x = {std::log(pv.Energy(0) * x_scaling), - std::log(pv.Energy(size - 1) * x_scaling)}; + grid.x = {std::log(bounds[0] * x_scaling), std::log(bounds[1] * x_scaling)}; grid.y.resize(size); - double delta - = fastpow(pv.Energy(size - 1) / pv.Energy(0), 1.0 / (size - 1)); + double delta_e = calc_log_delta(Span{bounds}); for (auto i : range(size)) { // Check that the grid has log spacing CELER_ASSERT(i == 0 - || soft_equal(delta, pv.Energy(i) / pv.Energy(i - 1))); + || soft_equal(delta_e, pv.Energy(i) / pv.Energy(i - 1))); grid.y[i] = pv[i] * y_scaling; } CELER_ENSURE(grid); diff --git a/src/corecel/cont/Span.hh b/src/corecel/cont/Span.hh index 8eb5f2c2ff..0e5aaa6f6c 100644 --- a/src/corecel/cont/Span.hh +++ b/src/corecel/cont/Span.hh @@ -188,8 +188,7 @@ class Span std::enable_if_t && (E2 == Extent || Extent == dynamic_extent), bool> = true> - CELER_CONSTEXPR_FUNCTION - Span(Span const& other) noexcept(ndebug_or_dyn) + CELER_CONSTEXPR_FUNCTION Span(Span other) noexcept(ndebug_or_dyn) : s_(other.data(), other.size()) { } @@ -204,7 +203,7 @@ class Span std::enable_if_t && Extent != dynamic_extent && E2 == dynamic_extent, bool> = true> - CELER_CONSTEXPR_FUNCTION explicit Span(Span const& other) noexcept( + CELER_CONSTEXPR_FUNCTION explicit Span(Span other) noexcept( ndebug_or_dyn) : s_(other.data(), other.size()) { diff --git a/src/corecel/grid/UniformGrid.hh b/src/corecel/grid/UniformGrid.hh index f3093c85ab..8197a05f5f 100644 --- a/src/corecel/grid/UniformGrid.hh +++ b/src/corecel/grid/UniformGrid.hh @@ -80,7 +80,7 @@ UniformGrid::UniformGrid(UniformGridData const& data) : data_(data) CELER_FUNCTION auto UniformGrid::operator[](size_type i) const -> value_type { CELER_EXPECT(i < data_.size); - return data_.front + data_.delta * i; + return data_.front + data_.delta * static_cast(i); } //---------------------------------------------------------------------------// diff --git a/src/corecel/grid/UniformGridData.hh b/src/corecel/grid/UniformGridData.hh index 0eb000d032..8a40a673a7 100644 --- a/src/corecel/grid/UniformGridData.hh +++ b/src/corecel/grid/UniformGridData.hh @@ -87,7 +87,8 @@ UniformGridData::from_bounds(EnumArray bounds, size_type size) result.size = size; result.front = bounds[Bound::lo]; result.back = bounds[Bound::hi]; - result.delta = (bounds[Bound::hi] - bounds[Bound::lo]) / (size - 1); + result.delta = (bounds[Bound::hi] - bounds[Bound::lo]) + / static_cast(size - 1); CELER_ENSURE(result); return result; } diff --git a/src/corecel/grid/VectorUtils.hh b/src/corecel/grid/VectorUtils.hh index cfd29d65fa..25321efa8e 100644 --- a/src/corecel/grid/VectorUtils.hh +++ b/src/corecel/grid/VectorUtils.hh @@ -68,7 +68,8 @@ template T calc_log_delta(Span grid) { CELER_EXPECT(grid.size() > 1); - return fastpow(grid.back() / grid.front(), T(1) / (grid.size() - 1)); + return fastpow(grid.back() / grid.front(), + T(1) / static_cast(grid.size() - 1)); } //---------------------------------------------------------------------------// diff --git a/src/geocel/rasterize/SafetyImager.t.hh b/src/geocel/rasterize/SafetyImager.t.hh index b673600007..d78625cf45 100644 --- a/src/geocel/rasterize/SafetyImager.t.hh +++ b/src/geocel/rasterize/SafetyImager.t.hh @@ -41,8 +41,9 @@ void SafetyImager::operator()(ImageParams const& image, std::string filename) out << nlohmann::json(image).dump() << std::endl; auto const& scalars = image.scalars(); - real_type max_distance = celeritas::max(scalars.dims[0], scalars.dims[1]) - * scalars.pixel_width; + auto max_distance = static_cast( + celeritas::max(scalars.dims[0], scalars.dims[1])) + * scalars.pixel_width; detail::SafetyCalculator calc_safety{ GeoTrackView{geo_->host_ref(), host_state_.ref(), TrackSlotId{0}}, diff --git a/test/corecel/random/Histogram.hh b/test/corecel/random/Histogram.hh index 523f970494..eed9858afb 100644 --- a/test/corecel/random/Histogram.hh +++ b/test/corecel/random/Histogram.hh @@ -103,7 +103,8 @@ void Histogram::operator()(double value) } else if (frac < 1.0) { - auto index = static_cast(frac * counts_.size()); + auto index = static_cast( + frac * static_cast(counts_.size())); CELER_ASSERT(index < counts_.size()); ++counts_[index]; } diff --git a/test/corecel/random/distribution/PoissonDistribution.test.cc b/test/corecel/random/distribution/PoissonDistribution.test.cc index 3764e85bff..633d2cc1d8 100644 --- a/test/corecel/random/distribution/PoissonDistribution.test.cc +++ b/test/corecel/random/distribution/PoissonDistribution.test.cc @@ -62,7 +62,8 @@ TEST(PoissonDistributionKnuthTest, large) HistogramSampler calc_histogram(8, {0, 50}, 1000); std::vector actual; // Test default lambda=1 - actual.push_back(calc_histogram(PoissonDistributionKnuth{})); + actual.push_back(calc_histogram(static_cast_double, + PoissonDistributionKnuth{})); for (auto i : range(1, 18)) { diff --git a/test/testdetail/TestMacrosImpl.cc b/test/testdetail/TestMacrosImpl.cc index 07d45c23a3..2fc60b27b8 100644 --- a/test/testdetail/TestMacrosImpl.cc +++ b/test/testdetail/TestMacrosImpl.cc @@ -50,12 +50,11 @@ int num_digits(unsigned long val) * * where too long means > digits digits. */ -char const* -trunc_string(unsigned int digits, char const* str, char const* trunc) +char const* trunc_string(int digits, char const* str, char const* trunc) { CELER_EXPECT(str && trunc); CELER_EXPECT(digits > 0); - CELER_EXPECT(std::strlen(trunc) <= digits); + CELER_EXPECT(std::strlen(trunc) <= static_cast(digits)); if (std::strlen(str) > digits) { diff --git a/test/testdetail/TestMacrosImpl.hh b/test/testdetail/TestMacrosImpl.hh index 8932ef4587..dc0c1f1226 100644 --- a/test/testdetail/TestMacrosImpl.hh +++ b/test/testdetail/TestMacrosImpl.hh @@ -31,8 +31,7 @@ namespace testdetail int num_digits(unsigned long val); // Return a replacement string if the given string is too long -char const* -trunc_string(unsigned int digits, char const* str, char const* trunc); +char const* trunc_string(int digits, char const* str, char const* trunc); //---------------------------------------------------------------------------// // INLINE DEFINITIONS @@ -367,8 +366,8 @@ IsRangeEqImpl(Iter1 e_iter, BinaryOp comp) { using size_type = std::size_t; - size_type expected_size = std::distance(e_iter, e_end); - size_type actual_size = std::distance(a_iter, a_end); + auto expected_size = std::distance(e_iter, e_end); + auto actual_size = std::distance(a_iter, a_end); // First, check that the sizes are equal if (expected_size != actual_size) From 4592d9a459fdbe6caaceb0054d3e220f167650aa Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Thu, 4 Jun 2026 12:27:59 -0400 Subject: [PATCH 25/28] fixup! Explicitly cast from int to double --- test/testdetail/TestMacrosImpl.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/testdetail/TestMacrosImpl.cc b/test/testdetail/TestMacrosImpl.cc index 2fc60b27b8..b2ee4dc796 100644 --- a/test/testdetail/TestMacrosImpl.cc +++ b/test/testdetail/TestMacrosImpl.cc @@ -56,7 +56,7 @@ char const* trunc_string(int digits, char const* str, char const* trunc) CELER_EXPECT(digits > 0); CELER_EXPECT(std::strlen(trunc) <= static_cast(digits)); - if (std::strlen(str) > digits) + if (std::strlen(str) > static_cast(digits)) { return trunc; } From c13e989d945936b8f117120942a9b316eb1d3f82 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Thu, 4 Jun 2026 16:13:54 -0400 Subject: [PATCH 26/28] fixup! Explicitly cast from int to double --- src/celeritas/ext/detail/GeantProcessImporter.cc | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/celeritas/ext/detail/GeantProcessImporter.cc b/src/celeritas/ext/detail/GeantProcessImporter.cc index 674c6a22f9..6cca34389a 100644 --- a/src/celeritas/ext/detail/GeantProcessImporter.cc +++ b/src/celeritas/ext/detail/GeantProcessImporter.cc @@ -35,9 +35,7 @@ #include "corecel/Assert.hh" #include "corecel/cont/Range.hh" -#include "corecel/cont/Span.hh" #include "corecel/data/HyperslabIndexer.hh" -#include "corecel/grid/VectorUtils.hh" #include "corecel/io/Logger.hh" #include "corecel/math/Algorithms.hh" #include "corecel/math/SoftEqual.hh" @@ -389,18 +387,20 @@ inp::UniformGrid import_physics_log_vector(G4PhysicsVector const& pv, double const x_scaling = native_value_from_clhep(units[0]); double const y_scaling = native_value_from_clhep(units[1]); auto size = pv.GetVectorLength(); - double const bounds[] = {pv.Energy(0), pv.Energy(size - 1)}; + CELER_ASSERT(size > 1); inp::UniformGrid grid; - grid.x = {std::log(bounds[0] * x_scaling), std::log(bounds[1] * x_scaling)}; + grid.x = {std::log(pv.Energy(0) * x_scaling), + std::log(pv.Energy(size - 1) * x_scaling)}; grid.y.resize(size); - double delta_e = calc_log_delta(Span{bounds}); + double delta = fastpow(pv.Energy(size - 1) / pv.Energy(0), + 1.0 / static_cast(size - 1)); for (auto i : range(size)) { // Check that the grid has log spacing CELER_ASSERT(i == 0 - || soft_equal(delta_e, pv.Energy(i) / pv.Energy(i - 1))); + || soft_equal(delta, pv.Energy(i) / pv.Energy(i - 1))); grid.y[i] = pv[i] * y_scaling; } CELER_ENSURE(grid); From 0a8dd50befbbdbf04c618b15b4d59d65d8fd022f Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Thu, 4 Jun 2026 16:55:26 -0400 Subject: [PATCH 27/28] Fix cuda warning --- src/corecel/random/distribution/TruncatedDistribution.hh | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/corecel/random/distribution/TruncatedDistribution.hh b/src/corecel/random/distribution/TruncatedDistribution.hh index 9fa85fad1b..3477fc8900 100644 --- a/src/corecel/random/distribution/TruncatedDistribution.hh +++ b/src/corecel/random/distribution/TruncatedDistribution.hh @@ -107,6 +107,11 @@ TruncatedDistribution::operator()(Generator& rng) -> result_type internal); } } + else + { + // CUDA 12.6 causes warning about unused variable + CELER_DISCARD(num_remaining_samples); + } result = sample_(rng); } while (result < lower_ || result > upper_); From 43089a2d3c5b338dcaf423306e77cd40c1b1c1db Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Thu, 4 Jun 2026 18:02:41 -0400 Subject: [PATCH 28/28] Fix abs/fabs --- src/celeritas/em/msc/detail/UrbanMscScatter.hh | 2 +- src/corecel/math/FerrariSolver.hh | 2 +- test/TestMacros.test.cc | 3 ++- test/celeritas/em/UrbanMsc.test.cc | 3 +-- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/celeritas/em/msc/detail/UrbanMscScatter.hh b/src/celeritas/em/msc/detail/UrbanMscScatter.hh index a57e887893..7d352201c0 100644 --- a/src/celeritas/em/msc/detail/UrbanMscScatter.hh +++ b/src/celeritas/em/msc/detail/UrbanMscScatter.hh @@ -526,7 +526,7 @@ UrbanMscScatter::compute_theta0(ParticleTrackView const& particle) const constexpr units::MevEnergy c_highland{13.6}; real_type theta0 = c_highland.value() - * std::abs(value_as(particle.charge())) + * std::fabs(value_as(particle.charge())) * std::sqrt(y) * invbetacp; // Correction factor from e- scattering data diff --git a/src/corecel/math/FerrariSolver.hh b/src/corecel/math/FerrariSolver.hh index 47a9291079..b5477e35d9 100644 --- a/src/corecel/math/FerrariSolver.hh +++ b/src/corecel/math/FerrariSolver.hh @@ -338,7 +338,7 @@ FerrariSolver::real_roots_normalized_cubic(real_type b, else // One real and two complex roots, solve for real root with Cardano { real_type nr_a = -signum(r) - * std::cbrt(std::abs(r) + std::sqrt(discrim)); + * std::cbrt(std::fabs(r) + std::sqrt(discrim)); real_type nr_b = nr_a == 0 ? 0 : q / nr_a; real_type z0 = nr_a + nr_b - third_b; return Real3(z0, no_solution_, no_solution_); diff --git a/test/TestMacros.test.cc b/test/TestMacros.test.cc index 7b022b97df..34ed70f0e2 100644 --- a/test/TestMacros.test.cc +++ b/test/TestMacros.test.cc @@ -6,6 +6,7 @@ //---------------------------------------------------------------------------// #include "TestMacros.hh" +#include #include #include "corecel/cont/Array.hh" @@ -445,7 +446,7 @@ inline ::testing::AssertionResult IsRefEq(char const* expr1, { ::celeritas::test::AssertionHelper result(expr1, expr2); - if (std::abs(val1.val - val2.val) >= tol.val) + if (std::fabs(val1.val - val2.val) >= tol.val) { result.fail() << " foo: " << val1.val << " != " << val2.val << " within " << tolexpr << " (" << tol.val << ")"; diff --git a/test/celeritas/em/UrbanMsc.test.cc b/test/celeritas/em/UrbanMsc.test.cc index 2332a00214..4904646814 100644 --- a/test/celeritas/em/UrbanMsc.test.cc +++ b/test/celeritas/em/UrbanMsc.test.cc @@ -297,8 +297,7 @@ TEST_F(UrbanMscTest, step_conversion) ASSERT_NO_THROW(true_step = geo_to_true(gp.step)); /* * TODO: large relative error -0.00081720192362734587 when - pstep - * is near or equal to range: + * pstep is near or equal to range: * z -> g: Low energy or range-limited step: slope = 1.6653345369377e-15