corecel/math/PdfUtils.hh

0001 //------------------------------- -*- C++ -*- -------------------------------//
0002 // Copyright Celeritas contributors: see top-level COPYRIGHT file for details
0003 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0004 //---------------------------------------------------------------------------//
0005 //! \file corecel/math/PdfUtils.hh
0006 //---------------------------------------------------------------------------//
0007 #pragma once
0008
0009 #include <utility>
0010
0011 #include "corecel/cont/Array.hh"
0012 #include "corecel/cont/Range.hh"
0013 #include "corecel/cont/Span.hh"
0014
0015 #include "Algorithms.hh"
0016
0017 namespace celeritas
0018 {
0019 //---------------------------------------------------------------------------//
0020 /*!
0021  * Calculate the integral of a piecewise rectangular function.
0022  *
0023  * The value at the left point is taken for the interval.
0024  *
0025  * \todo Remove rectangular integrator (currently and likely always unused) and
0026  * possibly add higher-order integration methods if needed (e.g. for cross
0027  * section generation)
0028  */
0029 struct PostRectangleSegmentIntegrator
0030 {
0031     template<class T>
0032     T operator()(Array<T, 2> lo, Array<T, 2> hi) const
0033     {
0034         return (hi[0] - lo[0]) * lo[1];
0035     }
0036 };
0037
0038 //---------------------------------------------------------------------------//
0039 /*!
0040  * Calculate the integral of a piecewise linear function.
0041  */
0042 struct TrapezoidSegmentIntegrator
0043 {
0044     template<class T>
0045     T operator()(Array<T, 2> lo, Array<T, 2> hi) const
0046     {
0047         return T(0.5) * (hi[0] - lo[0]) * (hi[1] + lo[1]);
0048     }
0049 };
0050
0051 //---------------------------------------------------------------------------//
0052 /*!
0053  * Integrate a piecewise function.
0054  *
0055  * To construct a CDF, `init` should be zero, and the destination should be
0056  * normalized by its final value afterward.
0057  */
0058 template<class I>
0059 class SegmentIntegrator
0060 {
0061   public:
0062     //! Construct with integrator
0063     explicit SegmentIntegrator(I&& integrate)
0064         : integrate_{std::forward<I>(integrate)}
0065     {
0066     }
0067
0068     //! Integrate a function
0069     template<class T, std::size_t N, std::size_t M>
0070     void operator()(Span<T const, N> x,
0071                     Span<T const, N> f,
0072                     Span<T, M> dst,
0073                     T init = {})
0074     {
0075         CELER_EXPECT(x.size() == f.size());
0076         CELER_EXPECT(x.size() == dst.size());
0077
0078         using Array2 = Array<T, 2>;
0079
0080         Array2 prev{x[0], f[0]};
0081         dst[0] = init;
0082         for (auto i : range(std::size_t{1}, x.size()))
0083         {
0084             Array2 cur{x[i], f[i]};
0085             init += integrate_(prev, cur);
0086             dst[i] = init;
0087             prev = cur;
0088         }
0089     }
0090
0091   private:
0092     I integrate_;
0093 };
0094
0095 //---------------------------------------------------------------------------//
0096 /*!
0097  * Estimate the mean and variance of a tabulated PDF.
0098  */
0099 class MomentCalculator
0100 {
0101   public:
0102     template<class T>
0103     struct Result
0104     {
0105         T mean{};
0106         T variance{};
0107     };
0108
0109   public:
0110     //! Estimate the mean and variance
0111     template<class T>
0112     Result<T> operator()(Span<T const> x, Span<T const> f)
0113     {
0114         CELER_EXPECT(x.size() == f.size());
0115         CELER_EXPECT(x.size() >= 2);
0116
0117         using Array2 = Array<T, 2>;
0118
0119         TrapezoidSegmentIntegrator integrate;
0120
0121         T integral{};
0122         T mean{};
0123         T variance{};
0124         Array2 prev{x[0], f[0]};
0125         for (auto i : range(std::size_t{1}, x.size()))
0126         {
0127             Array2 cur{x[i], f[i]};
0128             auto area = integrate(prev, cur);
0129             auto x_eval = T(0.5) * (cur[0] + prev[0]);
0130             integral += area;
0131             mean += area * x_eval;
0132             variance += area * ipow<2>(x_eval);
0133             prev = cur;
0134         }
0135         mean /= integral;
0136         variance = variance / integral - ipow<2>(mean);
0137         return {mean, variance};
0138     }
0139 };
0140
0141 //---------------------------------------------------------------------------//
0142 /*!
0143  * Normalize a vector by the final value and check for monotonicity.
0144  */
0145 template<class T, std::size_t N>
0146 inline void normalize_cdf(Span<T, N> x)
0147 {
0148     CELER_EXPECT(!x.empty());
0149     CELER_EXPECT(x.back() > 0);
0150     T norm{1 / x.back()};
0151     for (auto i : range(x.size() - 1))
0152     {
0153         CELER_ASSERT(x[i + 1] >= x[i]);
0154         x[i] *= norm;
0155     }
0156     x.back() = 1;
0157 }
0158
0159 //---------------------------------------------------------------------------//
0160 }  // namespace celeritas