Program Listing for File Extinct.hpp

↰ Return to documentation for file (src/extinct/Extinct.hpp)

#ifndef EXTINCT_HPP
#define EXTINCT_HPP

#include "../interpolation/Interpolator1D.hpp"
#include "../phot/FilterTabulated.hpp"
#include "../utils/Constants.hpp"
#include "../utils/HDF5Utils.hpp"
#include "../utils/TOMLUtils.hpp"
#include "hdf5.h" // NOLINT(misc-include-cleaner) -- see HDF5Utils.hpp's own comment on including hdf5.h wholesale
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <filesystem>
#include <stdexcept>
#include <string>
#include <toml.hpp>
#include <utility>
#include <vector>

namespace io
{
    class SimControls;
} // namespace io

namespace extinct
{
    inline static const std::string defaultRegistry = // NOLINT(bugprone-throwing-static-initialization,cert-err58-cpp) -- built from fixed string literals, so the (theoretically throwing) path conversion can never actually throw here
        (std::filesystem::path("data") / std::filesystem::path("extinct")
        / std::filesystem::path("extinct.toml"));
    inline static const std::string defaultVRegistry = // NOLINT(bugprone-throwing-static-initialization,cert-err58-cpp) -- see defaultRegistry above
        (std::filesystem::path("data") / std::filesystem::path("filters")
        / std::filesystem::path("V_filter.toml"));
    class Extinct
    {
    public:

        Extinct(const std::string& extinctName,
            const io::SimControls& controls,
            const std::string& registryName = defaultRegistry) :
            controls_(controls)
        {
            loadCurveImpl(extinctName, registryName);
        }

        // Copyable (rebinding controls_ to the same referent) but not
        // assignable (controls_ can't be reseated), matching Specsyn's
        // own identical copy/move declarations exactly -- see its
        // comment. Never actually copied/assigned in practice: every
        // Extinct is held via unique_ptr (see SimControls's own
        // extinct_).
        Extinct(const Extinct&) = default;
        Extinct(Extinct&&) = default;
        auto operator=(const Extinct&) -> Extinct& = delete;
        auto operator=(Extinct&&) -> Extinct& = delete;
        ~Extinct() = default;

        void loadCurve(const std::string& extinctName,
            const std::string& registryName = defaultRegistry)
        {
            Extinct tmp(*this);
            tmp.loadCurveImpl(extinctName, registryName);
            commitFrom(tmp);
        }

        void rebuildCache();

        // Observers

        [[nodiscard]] auto controls() const -> const io::SimControls& { return controls_; }

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

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

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

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

        [[nodiscard]] auto wlObs() const -> std::vector<double>;

        [[nodiscard]] auto wlOffset() const { return wlOffset_; }

        [[nodiscard]] auto applyExtinction(const double A_V, // NOLINT(readability-identifier-naming) -- see above
            const std::vector<double>& spec) const -> std::vector<double>
        {
            if (wl_.empty())
            {
                throw std::runtime_error(
                    "Extinct::applyExtinction: this Extinct has no wavelength "
                    "grid -- controls_.specsyn() was null the last time "
                    "rebuildCache() ran, so there is nothing to extinguish "
                    "spec against");
            }
            std::vector<double> result(wl_.size());
            for (std::size_t i = 0; i < wl_.size(); i++)
            {
                result[i] = spec[wlOffset_ + i] * std::exp(-A_V * extinct_[i]); // NOLINT(cppcoreguidelines-pro-bounds-avoid-unchecked-container-access) -- result and extinct_ both have size wl_.size(), and i is bounded by wl_.size(); spec has the same size as controls_.specsyn()->wl() (this class's own contract), so wlOffset_ + i stays in bounds
            }
            return result;
        }

        [[nodiscard]] auto applyExtinctionCts(const std::vector<double>& spec,
            const bool nebular = false) const -> std::vector<double>
        {
            if (wl_.empty())
            {
                throw std::runtime_error(
                    "Extinct::applyExtinctionCts: this Extinct has no "
                    "wavelength grid -- controls_.specsyn() was null the "
                    "last time rebuildCache() ran, so there is nothing to "
                    "extinguish spec against");
            }
            const auto& fac = nebular ? extinctionFacCtsNeb_ : extinctionFacCts_;
            std::vector<double> result(wl_.size());
            for (std::size_t i = 0; i < wl_.size(); i++)
            {
                result[i] = spec[wlOffset_ + i] * fac[i]; // NOLINT(cppcoreguidelines-pro-bounds-avoid-unchecked-container-access) -- result and fac (extinctionFacCts_ or extinctionFacCtsNeb_) all have size wl_.size() (by computeExtinctionFacCts()'s own contract), and i is bounded by wl_.size(); spec has the same size as the wl originally passed to the constructor, so wlOffset_ + i stays in bounds -- see applyExtinction()'s own identical indexing
            }
            return result;
        }

        [[nodiscard]] auto applyExtinctionLines(const double A_V, // NOLINT(readability-identifier-naming) -- see applyExtinction()'s own identical NOLINT
            const std::vector<double>& lineLum) const -> std::vector<double>
        {
            if (wl_.empty())
            {
                throw std::runtime_error(
                    "Extinct::applyExtinctionLines: this Extinct has no "
                    "wavelength grid -- controls_.specsyn() was null the "
                    "last time rebuildCache() ran, so extinctLines_ could "
                    "not be meaningfully normalized");
            }
            std::vector<double> result(extinctLines_.size());
            for (std::size_t i = 0; i < extinctLines_.size(); i++)
            {
                result[i] = lineLum[i] * std::exp(-A_V * extinctLines_[i]); // NOLINT(cppcoreguidelines-pro-bounds-avoid-unchecked-container-access) -- result and extinctLines_ both have size extinctLines_.size(), and i is bounded by extinctLines_.size(); lineLum has the same size, by this method's own contract
            }
            return result;
        }

        [[nodiscard]] auto applyExtinctionCtsLines(const std::vector<double>& lineLum) const -> std::vector<double>
        {
            if (wl_.empty())
            {
                throw std::runtime_error(
                    "Extinct::applyExtinctionCtsLines: this Extinct has no "
                    "wavelength grid -- controls_.specsyn() was null the "
                    "last time rebuildCache() ran, so extinctLines_ could "
                    "not be meaningfully normalized");
            }
            std::vector<double> result(extinctLines_.size());
            for (std::size_t i = 0; i < extinctLines_.size(); i++)
            {
                result[i] = lineLum[i] * extinctionFacCtsLines_[i]; // NOLINT(cppcoreguidelines-pro-bounds-avoid-unchecked-container-access) -- result and extinctionFacCtsLines_ both have size extinctLines_.size() (extinctionFacCtsLines_ by computeExtinctionFacCtsLines()'s own contract), and i is bounded by extinctLines_.size(); lineLum has the same size, by this method's own contract
            }
            return result;
        }

    private:

        void commitFrom(Extinct& tmp)
        {
            wlDat_ = std::move(tmp.wlDat_);
            extinctDat_ = std::move(tmp.extinctDat_);
            wl_ = std::move(tmp.wl_);
            extinct_ = std::move(tmp.extinct_);
            wlOffset_ = tmp.wlOffset_;
            extinctLines_ = std::move(tmp.extinctLines_);
            extinctionFacCts_ = std::move(tmp.extinctionFacCts_);
            extinctionFacCtsNeb_ = std::move(tmp.extinctionFacCtsNeb_);
            extinctionFacCtsLines_ = std::move(tmp.extinctionFacCtsLines_);
        }

        void loadCurveImpl(const std::string& extinctName,
            const std::string& registryName)
        {
            // Locate and parse the registry file
            const auto [registry, registryPath] =
                utils::parseTOMLFile(registryName, "Extinct");

            // Validate that the registry actually lists this curve
            const auto curves = utils::getStringArrayField(registry, "curves");
            if (std::ranges::find(curves, extinctName) == curves.end())
            {
                throw std::runtime_error(
                    "Extinct: registry " + registryPath.string() +
                    " has no extinction curve '" + extinctName + "'");
            }

            // The registry's top-level "file" entry names the HDF5 file
            // holding the actual curve data, relative to the directory
            // containing the registry itself
            const auto h5Name = registry["file"].value<std::string>(); // NOLINT(cppcoreguidelines-pro-bounds-avoid-unchecked-container-access) -- toml::table::operator[] is a keyed lookup, not a bounds-checkable container index; a missing key just yields a null node_view, handled by value<std::string>() returning nullopt
            if (!h5Name.has_value())
            {
                throw std::runtime_error(
                    "Extinct: registry " + registryPath.string() +
                    " is missing required 'file' field");
            }
            const auto h5Path = registryPath.parent_path() / h5Name.value();

            // NOLINTBEGIN(misc-include-cleaner) -- see HDF5Utils.hpp's own comment
            const hid_t file = H5Fopen(h5Path.string().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
            if (file < 0)
            {
                throw std::runtime_error(
                    "Extinct: unable to open HDF5 file " + h5Path.string());
            }
            const hid_t grp = H5Gopen2(file, extinctName.c_str(), H5P_DEFAULT);
            if (grp < 0)
            {
                H5Fclose(file);
                throw std::runtime_error(
                    "Extinct: unable to open group " + extinctName +
                    " in HDF5 file " + h5Path.string());
            }

            wlDat_ = utils::readDataset1D(grp, "wavelength", "Extinct");
            extinctDat_ = utils::readDataset1D(grp, "kappa", "Extinct");

            H5Gclose(grp);
            H5Fclose(file);
            // NOLINTEND(misc-include-cleaner)

            rebuildCacheImpl();
        }

        static void normalize(const std::vector<double>& wl,
            std::vector<double>& extinct,
            std::vector<double>& extinctLines,
            const std::string& vRegistry = defaultVRegistry)
        {
            const phot::FilterTabulated vFilt("Generic", "Johnson", "V", vRegistry);

            // V filter response, converted from wavelength to
            // frequency and reversed back into increasing order (wl
            // increases, so nu = c/wl decreases), then rebuilt as a
            // frequency-space interpolator
            std::vector<double> filtNu(vFilt.wl().size());
            std::ranges::transform(vFilt.wl(), filtNu.begin(),
                [](const double w) -> double { return utils::c / (w * utils::Angstrom); });
            std::ranges::reverse(filtNu);
            const std::vector<double> filtResp(vFilt.responseData().rbegin(), vFilt.responseData().rend());
            const interp::Interpolator1D<1> respInterp(filtNu, filtResp);

            // Same wavelength-to-frequency transformation for the
            // extinction curve
            std::vector<double> extNu(wl.size());
            std::ranges::transform(wl, extNu.begin(),
                [](const double w) -> double { return utils::c / (w * utils::Angstrom); });
            std::ranges::reverse(extNu);
            const std::vector<double> extRev(extinct.rbegin(), extinct.rend());
            const interp::Interpolator1D<1> extInterp(extNu, extRev);

            const double norm = (std::log(100.0) / 5.0) *
                respInterp.integ(respInterp.xMin(), respInterp.xMax()) /
                respInterp.integ(extInterp, respInterp.xMin(), respInterp.xMax());

            for (double& e : extinct) { e *= norm; }
            for (double& e : extinctLines) { e *= norm; }
        }

        [[nodiscard]] auto extinctFac(const double A_V) const -> std::vector<double> // NOLINT(readability-identifier-naming) -- see applyExtinction()'s own identical NOLINT
        {
            std::vector<double> result(wl_.size());
            for (std::size_t i = 0; i < wl_.size(); i++)
            {
                result[i] = std::exp(-A_V * extinct_[i]); // NOLINT(cppcoreguidelines-pro-bounds-avoid-unchecked-container-access) -- result and extinct_ both have size wl_.size(), and i is bounded by wl_.size()
            }
            return result;
        }

        [[nodiscard]] auto extinctFacLines(const double A_V) const -> std::vector<double> // NOLINT(readability-identifier-naming) -- see applyExtinction()'s own identical NOLINT
        {
            std::vector<double> result(extinctLines_.size());
            for (std::size_t i = 0; i < extinctLines_.size(); i++)
            {
                result[i] = std::exp(-A_V * extinctLines_[i]); // NOLINT(cppcoreguidelines-pro-bounds-avoid-unchecked-container-access) -- result and extinctLines_ both have size extinctLines_.size(), and i is bounded by extinctLines_.size()
            }
            return result;
        }

        void initExtinctLines(const interp::Interpolator1D<1>& interp);

        using ExtinctFacFn = std::vector<double> (Extinct::*)(double) const;

        [[nodiscard]] auto expectedExtinctFac(ExtinctFacFn fac, std::size_t nInt,
            bool nebular) const -> std::vector<double>;

        void computeExtinctionFacCts();

        void computeExtinctionFacCtsLines();

        void rebuildCacheImpl();

        const io::SimControls& controls_; // NOLINT(cppcoreguidelines-avoid-const-or-ref-data-members) -- deliberately a live reference, not a copy: see this member's own comment for why. Only ever used through the same non-copyable, non-movable ownership pattern (unique_ptr in SimControls's own extinct_) as every other class with a reference member in this codebase (e.g. Specsyn's own controls_), so the usual objection (disabling implicit copy/move assignment) doesn't apply in practice.

        std::vector<double> wlDat_;
        std::vector<double> extinctDat_;
        std::vector<double> wl_;
        std::vector<double> extinct_;
        std::size_t wlOffset_ = 0;
        std::vector<double> extinctLines_;

        std::vector<double> extinctionFacCts_;

        std::vector<double> extinctionFacCtsNeb_;

        std::vector<double> extinctionFacCtsLines_;
    };

} // namespace extinct

#endif // EXTINCT_HPP