Skip to content

Commit

Permalink
Add a helper class to find the grid value corresponding to a sampled …
Browse files Browse the repository at this point in the history
…CDF value
  • Loading branch information
amandalund committed Nov 28, 2024
1 parent 9fb13d4 commit 0d187ba
Show file tree
Hide file tree
Showing 4 changed files with 149 additions and 15 deletions.
84 changes: 84 additions & 0 deletions src/celeritas/grid/InverseCdfFinder.hh
Original file line number Diff line number Diff line change
@@ -0,0 +1,84 @@
//----------------------------------*-C++-*----------------------------------//
// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers.
// See the top-level COPYRIGHT file for details.
// SPDX-License-Identifier: (Apache-2.0 OR MIT)
//---------------------------------------------------------------------------//
//! \file celeritas/grid/InverseCdfFinder.hh
//---------------------------------------------------------------------------//
#pragma once

#include "corecel/Assert.hh"
#include "corecel/Macros.hh"
#include "corecel/cont/Range.hh"
#include "corecel/grid/Interpolator.hh"
#include "corecel/math/Algorithms.hh"

namespace celeritas
{
//---------------------------------------------------------------------------//
/*!
* Given a sampled CDF value, find the corresponding grid value.
*
* \tparam G Grid, e.g. \c UniformGrid or \c NonUniformGrid
* \tparam C Calculate the CDF at a given grid index
*
* Both the input grid and the CDF must be monotonically increasing. The
* sampled CDF value must be in range.
*/
template<class G, class C>
class InverseCdfFinder
{
public:
// Construct from grid and CDF calculator
inline CELER_FUNCTION InverseCdfFinder(G const& grid, C&& calc_cdf);

// Find and interpolate the grid value corresponding to the given CDF
inline CELER_FUNCTION real_type operator()(real_type cdf) const;

private:
G const& grid_;
C calc_cdf_;
};

//---------------------------------------------------------------------------//
// INLINE DEFINITIONS
//---------------------------------------------------------------------------//
/*!
* Construct from grid and CDF calculator.
*/
template<class G, class C>
CELER_FUNCTION
InverseCdfFinder<G, C>::InverseCdfFinder(G const& grid, C&& calc_cdf)
: grid_(grid), calc_cdf_(celeritas::forward<C>(calc_cdf))
{
}

//---------------------------------------------------------------------------//
/*!
* Find and interpolate the grid value corresponding to the given CDF.
*/
template<class G, class C>
CELER_FUNCTION real_type InverseCdfFinder<G, C>::operator()(real_type cdf) const
{
CELER_EXPECT(cdf >= 0 && cdf < 1);

// Find the grid index of the sampled CDF value
Range indices(grid_.size());
auto iter = celeritas::lower_bound(
indices.begin(), indices.end(), cdf, [this](size_type i, real_type c) {
return calc_cdf_(i) < c;
});
CELER_ASSERT(iter != indices.end());
if (calc_cdf_(*iter) != cdf)
{
--iter;
}
size_type i = iter - indices.begin();

// Calculate the grid value corresponding to the sampled CDF value
return LinearInterpolator<real_type>{
{calc_cdf_(i), grid_[i]}, {calc_cdf_(i + 1), grid_[i + 1]}}(cdf);
}

//---------------------------------------------------------------------------//
} // namespace celeritas
23 changes: 8 additions & 15 deletions src/celeritas/neutron/interactor/detail/CascadeCollider.hh
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@
#include "corecel/math/ArrayUtils.hh"
#include "celeritas/Quantities.hh"
#include "celeritas/Types.hh"
#include "celeritas/grid/InverseCdfFinder.hh"
#include "celeritas/phys/FourVector.hh"
#include "celeritas/random/distribution/GenerateCanonical.hh"
#include "celeritas/random/distribution/UniformRealDistribution.hh"
Expand Down Expand Up @@ -145,22 +146,14 @@ CELER_FUNCTION auto CascadeCollider::operator()(Engine& rng) -> FinalState

if (kin_energy_ < energy_grid.back())
{
// Find cos\theta from tabulated angular data for a given c.d.f.
// Find cos\theta from tabulated angular data for a given CDF
Grid cos_grid(cdf_grid.y, shared_.reals);
TwodGridCalculator calc_cdf(cdf_grid, shared_.reals);

size_type idx = cos_grid.size() - 2;
real_type cdf_upper = 0;
real_type cdf_lower = 1;

do
{
cdf_upper = cdf_lower;
cdf_lower = calc_cdf({kin_energy_, cos_grid[idx]});
} while (cdf_lower > cdf && idx-- > 0);

real_type frac = (cdf - cdf_lower) / (cdf_upper - cdf_lower);
cos_theta = fma(frac, cos_grid[idx + 1] - cos_grid[idx], cos_grid[idx]);
auto calc_cdf
= TwodGridCalculator(cdf_grid, shared_.reals)(kin_energy_);
cos_theta
= InverseCdfFinder(cos_grid, [&calc_cdf, &cos_grid](size_type i) {
return calc_cdf(cos_grid[i]);
})(cdf);
}
else
{
Expand Down
1 change: 1 addition & 0 deletions test/celeritas/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -286,6 +286,7 @@ celeritas_add_test(grid/GenericCalculator.test.cc)
celeritas_add_test(grid/GenericGridBuilder.test.cc)
celeritas_add_test(grid/GenericGridInserter.test.cc)
celeritas_add_test(grid/GridIdFinder.test.cc)
celeritas_add_test(grid/InverseCdfFinder.test.cc)
celeritas_add_test(grid/InverseRangeCalculator.test.cc)
celeritas_add_test(grid/PolyEvaluator.test.cc)
celeritas_add_test(grid/RangeCalculator.test.cc)
Expand Down
56 changes: 56 additions & 0 deletions test/celeritas/grid/InverseCdfFinder.test.cc
Original file line number Diff line number Diff line change
@@ -0,0 +1,56 @@
//----------------------------------*-C++-*----------------------------------//
// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers.
// See the top-level COPYRIGHT file for details.
// SPDX-License-Identifier: (Apache-2.0 OR MIT)
//---------------------------------------------------------------------------//
//! \file celeritas/grid/InverseCdfFinder.test.cc
//---------------------------------------------------------------------------//
#include "celeritas/grid/InverseCdfFinder.hh"

#include "corecel/data/Collection.hh"
#include "corecel/data/CollectionBuilder.hh"
#include "corecel/grid/NonuniformGrid.hh"

#include "celeritas_test.hh"

namespace celeritas
{
namespace test
{
//---------------------------------------------------------------------------//

class InverseCdfFinderTest : public Test
{
protected:
using GridT = NonuniformGrid<real_type>;

void SetUp() override
{
CollectionBuilder build(&data);
values = build.insert_back({0.0, 0.2, 0.5, 1.5, 6.0, 10.0});
ref = data;
}

ItemRange<real_type> values;
Collection<real_type, Ownership::value, MemSpace::host> data;
Collection<real_type, Ownership::const_reference, MemSpace::host> ref;
};

TEST_F(InverseCdfFinderTest, all)
{
GridT grid(values, ref);
std::vector<real_type> cdf{0.0, 0.1, 0.2, 0.8, 0.9, 1.0};
InverseCdfFinder find(grid, [&cdf](size_type i) { return cdf[i]; });

EXPECT_SOFT_EQ(0, find(0));
EXPECT_SOFT_EQ(0.1, find(0.05));
EXPECT_SOFT_EQ(0.2, find(0.1));
EXPECT_SOFT_EQ(1, find(0.5));
EXPECT_SOFT_EQ(6.0, find(0.9));
EXPECT_SOFT_EQ(8.0, find(0.95));
EXPECT_SOFT_EQ(9.99996, find(0.999999));
}

//---------------------------------------------------------------------------//
} // namespace test
} // namespace celeritas

0 comments on commit 0d187ba

Please sign in to comment.