Program Listing for File YieldChannel.hpp

↰ Return to documentation for file (src/yields/YieldChannel.hpp)

#ifndef YIELDCHANNEL_HPP
#define YIELDCHANNEL_HPP

#include "../elem/ElemCommons.hpp"
#include "../utils/GridBracket.hpp"
#include "../utils/ThreadVec.hpp"
#include "YieldCommons.hpp"
#include <algorithm>
#include <cassert>
#include <cstddef>
#include <cstdint>
#include <mdspan> // NOLINT(misc-include-cleaner)
#include <optional>
#include <stdexcept>
#include <string>
#include <vector>

namespace yields
{
    class YieldChannel
    {
    public:

        using Array3D = std::mdspan<const double, std::dextents<std::size_t, 3>>; // NOLINT(misc-include-cleaner)

        using ActiveArray2D = std::mdspan<const std::uint8_t, std::dextents<std::size_t, 2>>; // NOLINT(misc-include-cleaner)

        // Constructors and destructors

        YieldChannel(
            const YieldChannelDescriptor& descriptor,
            double fehMin,
            double fehMax,
            const std::string& registryName = defaultRegistry);

        ~YieldChannel() = default;

        // Neither copyable nor movable: massCache_/fehCache_ are
        // utils::ThreadVec, itself neither copyable nor movable (see
        // its own comment), for the same reason SpecsynLib's own,
        // identical dim1Cache_/dim2Cache_/dim3Cache_ leave SpecsynLib
        // neither copyable nor movable either -- a per-thread cache
        // has no well-defined "the copy/move target's own cache"
        // to become. Constructed once (an HDF5 read) and queried
        // through a reference/pointer thereafter, the same as
        // SpecsynLib itself.
        YieldChannel(const YieldChannel&) = delete;
        YieldChannel(YieldChannel&&) = delete;
        auto operator=(const YieldChannel&) -> YieldChannel& = delete;
        auto operator=(YieldChannel&&) -> YieldChannel& = delete;

        void rebuildYieldGrid(
            std::optional<double> mMin = std::nullopt,
            std::optional<double> mMax = std::nullopt,
            elem::IsotopeList isotopes = {},
            const elem::IsotopeList& skippedIsotopes = {});

        // Observers

        [[nodiscard]] auto channel() const { return descriptor_.channel_; }

        [[nodiscard]] auto descriptor() const -> YieldChannelDescriptor { return descriptor_; }

        [[nodiscard]] auto masses() const -> const std::vector<double>& { return masses_; }

        [[nodiscard]] auto massesOrig() const -> const std::vector<double>& { return massesOrig_; }

        [[nodiscard]] auto isotopes() const -> const elem::IsotopeList& { return isotopes_; }

        [[nodiscard]] auto isotopesOrig() const -> const elem::IsotopeList& { return isotopesOrig_; }

        [[nodiscard]] auto feH() const -> const std::vector<double>& { return feH_; }

        [[nodiscard]] auto requestedFehMin() const -> double { return requestedFehMin_; }

        [[nodiscard]] auto requestedFehMax() const -> double { return requestedFehMax_; }

        [[nodiscard]] auto yld() const -> Array3D
        {
            return Array3D(yieldData_.data(), feH_.size(), masses_.size(), isotopes_.size());
        }

        [[nodiscard]] auto hasYield(const double mass) const -> bool
        {
            return !masses_.empty() && mass >= masses_.front() && mass <= masses_.back();
        }

        [[nodiscard]] auto hasYield(const double mass, const double feH) const -> bool
        {
            if (!hasYield(mass)) { return false; }
            const double feHClamped = std::clamp(feH, feH_.front(), feH_.back());
            const auto bm = utils::findBracket(masses_, mass, massCache_());
            const auto bf = utils::findBracket(feH_, feHClamped, fehCache_());
            const ActiveArray2D active(yieldActive_.data(), feH_.size(), masses_.size());
            // NOLINTBEGIN(cppcoreguidelines-pro-bounds-avoid-unchecked-container-access) -- bm/bf indices are all < the corresponding grid's size by construction
            return active[bf.lo_, bm.lo_] != 0 || active[bf.lo_, bm.hi_] != 0 ||
                active[bf.hi_, bm.lo_] != 0 || active[bf.hi_, bm.hi_] != 0;
            // NOLINTEND(cppcoreguidelines-pro-bounds-avoid-unchecked-container-access)
        }

        [[nodiscard]] auto yield(const double mass, const double feH) const -> std::vector<double>
        {
            if (yieldData_.empty())
            {
                throw std::runtime_error(
                    "YieldChannel::yield: yieldData_ is empty -- call "
                    "rebuildYieldGrid() at least once before calling yield()");
            }
            assert(hasYield(mass));
            assert(feH >= feH_.front() && feH <= feH_.back());

            const auto bm = utils::findBracket(masses_, mass, massCache_());
            const auto bf = utils::findBracket(feH_, feH, fehCache_());
            const auto view = yld();
            const std::size_t niso = isotopes_.size();

            std::vector<double> result(niso, 0.0);
            // NOLINTBEGIN(cppcoreguidelines-pro-bounds-avoid-unchecked-container-access) -- bm/bf indices are all < the corresponding grid's size by construction, and this loop is a hot path where the cost of bounds checking matters -- see this method's own comment
            for (int fi = 0; fi < 2; ++fi)
            {
                const std::size_t f = (fi == 0) ? bf.lo_ : bf.hi_;
                const double wf = (fi == 0) ? (1.0 - bf.t_) : bf.t_;
                if (wf == 0.0) { continue; } // degenerate/singular feH axis or exact grid hit
                for (int mi = 0; mi < 2; ++mi)
                {
                    const std::size_t m = (mi == 0) ? bm.lo_ : bm.hi_;
                    const double wm = (mi == 0) ? (1.0 - bm.t_) : bm.t_;
                    const double weight = wf * wm;
                    if (weight == 0.0) { continue; } // exact grid hit on mass
                    for (std::size_t iso = 0; iso < niso; ++iso)
                    {
                        result[iso] += weight * view[f, m, iso];
                    }
                }
            }
            // NOLINTEND(cppcoreguidelines-pro-bounds-avoid-unchecked-container-access)
            return result;
        }

    private:

        double requestedFehMin_;
        double requestedFehMax_;
        YieldChannelDescriptor descriptor_;
        std::vector<double> masses_;
        std::vector<double> massesOrig_;
        elem::IsotopeList isotopes_;
        elem::IsotopeList isotopesOrig_;
        std::vector<double> feH_;           // NOLINT(readability-identifier-naming)
        std::vector<double> yieldData_;
        std::vector<double> yieldDataOrig_;
        std::vector<std::uint8_t> yieldActiveOrig_;

        std::vector<std::uint8_t> yieldActive_;
        mutable utils::ThreadVec<std::size_t> massCache_;
        mutable utils::ThreadVec<std::size_t> fehCache_;
    };

} // namespace yields

#endif // YIELDCHANNEL_HPP