Program Listing for File SimControls.hpp

↰ Return to documentation for file (src/io/SimControls.hpp)

#ifndef SIMCONTROLS_HPP
#define SIMCONTROLS_HPP

#include "../extinct/Extinct.hpp"
#include "../feedback/FeedbackCommons.hpp"
#include "../feedback/Winds.hpp"
#include "../nebular/Nebular.hpp"
#include "../nebular/NebularCommons.hpp"
#include "../pdfs/PDF.hpp"
#include "../phot/FilterCollection.hpp"
#include "../specsyn/Specsyn.hpp"
#include "../tracks/Tracks2D.hpp"
#include "../tracks/Tracks3D.hpp"
#include "../utils/ParseUtils.hpp"
#include "../utils/TrackedDeck.hpp"
#include "../yields/YieldCommons.hpp"
#include "../yields/Yields.hpp"
#include <cstddef>
#include <cstdint>
#include <filesystem>
#include <memory>
#include <stdexcept>
#include <string>
#include <toml.hpp>
#include <utility>
#include <vector>

namespace io
{
    inline const std::string imfPrefix = // 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("imfs")).string();

    class SimControls
    {
    public:

        enum class OutputMode : std::uint8_t {
            h5,
            h5divided,
            ascii
        };

        enum class SimType : std::uint8_t {
            cluster,
            galaxy,
            none
        };

        SimControls() = default;

        explicit SimControls(const toml::table& inputDeck);

        // Non-copyable and non-movable: the spectral synthesizer
        // readSpectra() builds, and the nebular emission grid
        // readNebular() builds, each store a live reference back to
        // *this (see Specsyn's own controls_ member, and Nebular's own
        // simControls_ member) -- this object's address must never
        // change after construction for those references to stay
        // valid. (specsyn_/filters_/extinct_/nebular_/yields_/tracks_/
        // constFeHTracks_ are all held by shared_ptr, so none of them
        // force non-copyability/non-movability on their own the way an
        // owned non-copyable value member would -- the live
        // back-reference above is the only reason this class can't be
        // either.)
        SimControls(const SimControls&) = delete;
        auto operator=(const SimControls&) -> SimControls& = delete;
        SimControls(SimControls&&) = delete;
        auto operator=(SimControls&&) -> SimControls& = delete;
        ~SimControls() = default;

        // Getters for control flow

        [[nodiscard]] auto simType() const { return simType_; }

        [[nodiscard]] auto outputMode() const { return outputMode_; }

        [[nodiscard]] auto modelName() const { return modelName_; }

        [[nodiscard]] auto outDir() const { return outDir_; }

        [[nodiscard]] auto verbosity() const { return verbosity_; }

        [[nodiscard]] auto strictInput() const { return strictInput_; }

        [[nodiscard]] auto unusedKeys() const -> const std::vector<utils::DeckKeyReport>& { return unusedKeys_; }

        [[nodiscard]] auto ignoredKeys() const -> const std::vector<utils::DeckKeyReport>& { return ignoredKeys_; }

        [[nodiscard]] auto nTrial() const { return nTrial_; }

        [[nodiscard]] auto outTimes() const -> std::vector<double>
        {
            if (!outTimes_.empty()) { return outTimes_; }
            return { outTimeDist_.draw() };
        }

        [[nodiscard]] auto intRelTol() const { return intRelTol_; }

        [[nodiscard]] auto intAbsTol() const { return intAbsTol_; }

        [[nodiscard]] auto intMaxIter() const { return intMaxIter_; }

        [[nodiscard]] auto z() const { return z_; }

        [[nodiscard]] auto checkpointInterval() const { return checkpointInterval_; }

        [[nodiscard]] auto inputDeckStr() const -> const std::string& { return inputDeckStr_; }

        [[nodiscard]] auto writeCluster() const { return writeCluster_; }

        [[nodiscard]] auto writeClusterSpec() const { return writeClusterSpec_; }

        [[nodiscard]] auto writeClusterPhot() const { return writeClusterPhot_; }

        [[nodiscard]] auto writeGalaxy() const { return writeGalaxy_; }

        [[nodiscard]] auto writeGalaxySpec() const { return writeGalaxySpec_; }

        [[nodiscard]] auto writeGalaxyPhot() const { return writeGalaxyPhot_; }

        [[nodiscard]] auto writeClusterYields() const { return writeClusterYields_; }

        [[nodiscard]] auto writeGalaxyYields() const { return writeGalaxyYields_; }

        [[nodiscard]] auto writeClusterFeedback() const { return writeClusterFeedback_; }

        [[nodiscard]] auto writeGalaxyFeedback() const { return writeGalaxyFeedback_; }

        void setIntRelTol(double tol) { intRelTol_ = tol; }

        void setIntAbsTol(double tol) { intAbsTol_ = tol; }

        void setIntMaxIter(std::size_t n) { intMaxIter_ = n; }

        void setZ(double z) { z_ = z; }

        void setCheckpointInterval(unsigned long interval) { checkpointInterval_ = interval; }

        void setModelName(std::string name) { modelName_ = std::move(name); }

        void setOutDir(std::string dir) { outDir_ = std::move(dir); }

        void setNTrial(unsigned long n) { nTrial_ = n; }

        void setOutputMode(OutputMode mode) { outputMode_ = mode; }

        void setVerbosity(unsigned int verbosity) { verbosity_ = verbosity; }

        void setOutTimes(std::vector<double> times)
        {
            if (times.empty())
            {
                throw std::invalid_argument(
                    "setOutTimes: output time list must not be empty");
            }
            outTimes_ = std::move(times);
            outTimeDist_ = pdfs::PDF();
        }

        void setWriteCluster(bool value) { writeCluster_ = value; }

        void setWriteClusterSpec(bool value) { writeClusterSpec_ = value; }

        void setWriteClusterPhot(bool value) { writeClusterPhot_ = value; }

        void setWriteGalaxy(bool value) { writeGalaxy_ = value; }

        void setWriteGalaxySpec(bool value) { writeGalaxySpec_ = value; }

        void setWriteGalaxyPhot(bool value) { writeGalaxyPhot_ = value; }

        void setWriteClusterYields(bool value) { writeClusterYields_ = value; }

        void setWriteGalaxyYields(bool value) { writeGalaxyYields_ = value; }

        void setWriteClusterFeedback(bool value) { writeClusterFeedback_ = value; }

        void setWriteGalaxyFeedback(bool value) { writeGalaxyFeedback_ = value; }

        void setYieldsChannelDecomposed(bool value) { yieldsChannelDecomposed_ = value; }

        void setNoDecay(bool value) { noDecay_ = value; }

        void setMinIsotopeLifetime(double value);

        // Getters for the physics settings
        [[nodiscard]] auto imf() const -> const auto& { return imf_; }

        [[nodiscard]] auto cmf() const -> const auto& { return cmf_; }

        [[nodiscard]] auto fehDist() const -> const auto& { return fehDist_; }

        [[nodiscard]] auto constFeH() const { return fehDist_.getMin() == fehDist_.getMax(); }

        [[nodiscard]] auto sfr() const -> const auto& { return sfr_; }

        [[nodiscard]] auto sfrDist() const -> const auto& { return sfrDist_; }

        [[nodiscard]] auto clf() const -> const auto& { return clf_; }

        [[nodiscard]] auto fCluster() const { return fCluster_; }

        [[nodiscard]] auto tracks() const -> std::shared_ptr<const tracks::Tracks3D> { return tracks_; }

        [[nodiscard]] auto tracks2D() const -> std::shared_ptr<const tracks::Tracks2D> { return constFeHTracks_; }

        [[nodiscard]] auto minStochMass() const { return minStochMass_; }

        [[nodiscard]] auto fracStochMass() const { return fracStochMass_; }

        [[nodiscard]] auto nonStochIMFMass() const -> double
        {
            const double mMin = imf_.getMin();
            if (minStochMass_ <= mMin) { return 0.0; }
            // If the whole IMF is non-stochastic, use its whole-range
            // integral and mean directly: PDF::integral(a, b) is 0 for
            // any zero-width [a, b], which would otherwise give 0 for a
            // delta-function IMF (zero-width support, but all of its
            // stars in it)
            if (minStochMass_ >= imf_.getMax()) { return imf_.expectationValue() * imf_.integral(); }
            return imf_.expectationValue(mMin, minStochMass_) * imf_.integral(mMin, minStochMass_);
        }

        [[nodiscard]] auto specsyn() const -> std::shared_ptr<const specsyn::Specsyn> { return specsyn_; }

        [[nodiscard]] auto filters() const -> std::shared_ptr<const phot::FilterCollection> { return filters_; }

        [[nodiscard]] auto avDist() const -> const auto& { return avDist_; }

        [[nodiscard]] auto avDistField() const -> const auto& { return avDistField_; }

        [[nodiscard]] auto avNebFac() const -> const auto& { return avNebFac_; }

        [[nodiscard]] auto outTimeDist() const -> const auto& { return outTimeDist_; }

        [[nodiscard]] auto allPDFs() const -> std::vector<std::pair<std::string, const pdfs::PDF*>>
        {
            return {
                {"imf", &imf_},
                {"cmf", &cmf_},
                {"fehDist", &fehDist_},
                {"sfr", &sfr_},
                {"sfrDist", &sfrDist_},
                {"clf", &clf_},
                {"avDist", &avDist_},
                {"avDistField", &avDistField_},
                {"avNebFac", &avNebFac_},
                {"outTimeDist", &outTimeDist_},
            };
        }

        [[nodiscard]] auto extinct() const -> std::shared_ptr<const extinct::Extinct> { return extinct_; }

        [[nodiscard]] auto yieldChannels() const -> const std::vector<yields::YieldChannelDescriptor>&
        {
            return yieldChannels_;
        }

        [[nodiscard]] auto yields() const -> std::shared_ptr<const yields::Yields> { return yields_; }

        [[nodiscard]] auto yieldsChannelDecomposed() const { return yieldsChannelDecomposed_; }

        [[nodiscard]] auto noDecay() const { return noDecay_; }

        [[nodiscard]] auto minIsotopeLifetime() const { return minIsotopeLifetime_; }

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

        void setSNMassLimits(std::vector<double> limits);

        [[nodiscard]] auto hasSN(double mass) const -> bool;

        [[nodiscard]] auto hasSN(double mass, double feH) const -> bool;

        [[nodiscard]] auto wrWindModel() const -> feedback::WRwindModel { return wrWindModel_; }

        void setWRWindModel(const feedback::WRwindModel model) { wrWindModel_ = model; }

        [[nodiscard]] auto obWindModel() const -> feedback::OBwindModel { return obWindModel_; }

        void setOBWindModel(const feedback::OBwindModel model) { obWindModel_ = model; }

        [[nodiscard]] auto agbWindModel() const -> feedback::AGBwindModel { return agbWindModel_; }

        void setAGBWindModel(const feedback::AGBwindModel model) { agbWindModel_ = model; }

        [[nodiscard]] auto otherWindModel() const -> feedback::OtherwindModel { return otherWindModel_; }

        void setOtherWindModel(const feedback::OtherwindModel model) { otherWindModel_ = model; }

        [[nodiscard]] auto winds() const -> std::shared_ptr<const feedback::Winds> { return winds_; }

        [[nodiscard]] auto nebControls() const -> const nebular::NebularControls& { return nebControls_; }

        [[nodiscard]] auto nebular() const -> std::shared_ptr<const nebular::Nebular> { return nebular_; }

        [[nodiscard]] auto computeLbol() const { return computeLbol_; }

        void setComputeLbol(bool value) { computeLbol_ = value; }

        // Physics-setting setters: thin wrappers around
        // utils::initPDFFromString(), letting a caller (e.g. from
        // Python, building a SimControls without an input deck at all)
        // set each physics PDF directly from a string playing the same
        // role as the corresponding toml key's value would in the
        // constructor -- a numerical value becomes a delta function,
        // anything else is interpreted as a PDF file name. See
        // utils::initPDFFromString()'s own comment for the exact rules.

        void setIMF(const std::string& imf);

        void setCMF(const std::string& cmf) { cmf_ = utils::initPDFFromString(cmf); }

        void setFeH(const std::string& feH);

        void setAVDist(const std::string& avDist);

        void setAVDistField(const std::string& avDistField);

        void setAVNebFac(const std::string& avNebFac);

        void setCLF(const std::string& clf) { clf_ = utils::initPDFFromString(clf); }

        void setSFR(const std::string& sfr);

        void setSFRDist(const std::string& sfrDist)
        {
            sfrDist_ = utils::initPDFFromString(sfrDist);
            sfr_ = pdfs::PDF();
        }

        [[nodiscard]] static auto buildConstantSFR(double sfr) -> pdfs::PDF;

        void setFCluster(double fCluster) { fCluster_ = fCluster; }

        void setNebControls(const nebular::NebularControls& nebControls)
        {
            auto updated = nebControls;
            updated.computeNeb_ = nebControls_.computeNeb_;
            nebControls_ = updated;
        }

        // Object-replacement setters: unlike the string-driven setters
        // above, these accept an already-built object -- e.g. from
        // Python, where Specsyn, FilterCollection, Tracks3D, Extinct,
        // and Nebular are all directly constructible -- and install it
        // in place of whatever this SimControls already holds. Of
        // these, Specsyn/Extinct/Nebular each store a live reference
        // to whichever SimControls they were built against (see each
        // class's own controls_/simControls_ comment), so their own
        // setters additionally validate that the object being
        // installed was actually built against *this, throwing
        // std::invalid_argument otherwise -- FilterCollection/
        // Tracks3D hold no such reference, so setFilters()/
        // setTracks() have nothing to validate.

        void setSpecsyn(std::unique_ptr<specsyn::Specsyn> specsyn)
        {
            if (specsyn) { checkSpecsynReplacement(*specsyn); }

            if (!extinct_)
            {
                specsyn_ = std::move(specsyn);
                return;
            }

            auto oldSpecsyn = std::move(specsyn_);
            specsyn_ = std::move(specsyn);
            try
            {
                extinct_->rebuildCache();
            }
            catch (...)
            {
                specsyn_ = std::move(oldSpecsyn);
                throw;
            }
        }

        void setFilters(std::unique_ptr<phot::FilterCollection> filters) { filters_ = std::move(filters); }

        void setTracks(std::unique_ptr<tracks::Tracks3D> tracks);

        void checkTracksCoverFeH(const tracks::Tracks3D& tracks) const;

        void checkTracksCoverIMF(const tracks::Tracks3D& tracks) const;

        void setExtinct(std::unique_ptr<extinct::Extinct> extinct)
        {
            if (extinct && &extinct->controls() != this)
            {
                throw std::invalid_argument(
                    "SimControls::setExtinct: extinct was constructed "
                    "against a different SimControls than this one -- "
                    "construct it with this same SimControls instead");
            }
            extinct_ = std::move(extinct);
        }

        void setNebular(std::unique_ptr<nebular::Nebular> nebular)
        {
            if (nebular && &nebular->simControls() != this)
            {
                throw std::invalid_argument(
                    "SimControls::setNebular: nebular was constructed "
                    "against a different SimControls than this one -- "
                    "construct it with this same SimControls instead");
            }
            nebular_ = std::move(nebular);
        }

        void setYields(std::unique_ptr<yields::Yields> yields)
        {
            if (yields) { checkYieldsReplacement(*yields); }
            yields_ = std::move(yields);
        }

        void setWinds(std::unique_ptr<feedback::Winds> winds)
        {
            if (winds) { checkWindsReplacement(*winds); }
            winds_ = std::move(winds);
        }

        void checkSpecsynReplacement(const specsyn::Specsyn& specsyn) const
        {
            if (&specsyn.controls() != this)
            {
                throw std::invalid_argument(
                    "SimControls::setSpecsyn: specsyn was constructed "
                    "against a different SimControls than this one -- "
                    "construct it with this same SimControls instead");
            }
            if (specsyn_)
            {
                checkFehRangeContains("SimControls::setSpecsyn", "spectral synthesizer",
                    specsyn.requestedFehMin(), specsyn.requestedFehMax(),
                    specsyn_->requestedFehMin(), specsyn_->requestedFehMax(),
                    "the current spectral synthesizer's own");
            }
            else if (fehDist_.valid())
            {
                checkFehRangeContains("SimControls::setSpecsyn", "spectral synthesizer",
                    specsyn.requestedFehMin(), specsyn.requestedFehMax(),
                    fehDist_.getMin(), fehDist_.getMax(), "the current [Fe/H] distribution");
            }
        }

        void checkYieldsReplacement(const yields::Yields& yields) const
        {
            if (&yields.controls() != this)
            {
                throw std::invalid_argument(
                    "SimControls::setYields: yields was constructed "
                    "against a different SimControls than this one -- "
                    "construct it with this same SimControls instead");
            }
            if (yields_)
            {
                checkFehRangeContains("SimControls::setYields", "Yields",
                    yields.requestedFehMin(), yields.requestedFehMax(),
                    yields_->requestedFehMin(), yields_->requestedFehMax(),
                    "the current Yields' own");
            }
            else if (fehDist_.valid())
            {
                checkFehRangeContains("SimControls::setYields", "Yields",
                    yields.requestedFehMin(), yields.requestedFehMax(),
                    fehDist_.getMin(), fehDist_.getMax(), "the current [Fe/H] distribution");
            }
        }

        void checkWindsReplacement(const feedback::Winds& winds) const
        {
            if (&winds.controls() != this)
            {
                throw std::invalid_argument(
                    "SimControls::setWinds: winds was constructed "
                    "against a different SimControls than this one -- "
                    "construct it with this same SimControls instead");
            }
        }

        void setMinStochMass(double minStochMass)
        {
            minStochMass_ = minStochMass;
            updateFracStochMass();
        }

    private:

        static void checkFehRangeContains(const std::string& who, const std::string& what,
            const double newMin, const double newMax,
            const double refMin, const double refMax, const std::string& refWhat)
        {
            if (!(newMin <= refMin && refMax <= newMax)) // NOLINT(readability-simplify-boolean-expr) -- the De Morgan form would accept a NaN range, since every comparison with NaN is false
            {
                throw std::invalid_argument(
                    who + ": the new " + what + " was constructed for [Fe/H] in [" +
                    std::to_string(newMin) + ", " + std::to_string(newMax) + "], which "
                    "does not contain " + refWhat + ", [" + std::to_string(refMin) + ", " +
                    std::to_string(refMax) + "]. Construct it for a range covering that.");
            }
        }

        static void checkIMFWithinTracks(const std::string& who, const pdfs::PDF& imf,
            const tracks::Tracks3D& tracks);

        void updateFracStochMass()
        {
            const double mMin = imf_.getMin();
            const double mMax = imf_.getMax();
            if (minStochMass_ <= mMin) { fracStochMass_ = 1.0; return; }
            if (minStochMass_ >= mMax) { fracStochMass_ = 0.0; return; }
            const double massTot = imf_.expectationValue(mMin, mMax);
            const double massStoch = imf_.expectationValue(minStochMass_, mMax) *
                imf_.integral(minStochMass_, mMax);
            fracStochMass_ = massStoch / massTot;
        }

        void initControlFlow(const utils::TrackedDeck& inputDeck);

        void initPhysics(const utils::TrackedDeck& inputDeck);

        void setOutputTimes(const utils::TrackedDeck& inputDeck);

        void readOutput(const utils::TrackedDeck& inputDeck);

        void readTracks(const utils::TrackedDeck& inputDeck);

        void readSpectra(const utils::TrackedDeck& inputDeck);

        void readFilters(const utils::TrackedDeck& inputDeck);

        void readExtinct(const utils::TrackedDeck& inputDeck);

        void readYields(const utils::TrackedDeck& inputDeck);

        void readFeedback(const utils::TrackedDeck& inputDeck);

        void readNebular(const utils::TrackedDeck& inputDeck);

        void reportUnusedKeys(const utils::TrackedDeck& deck);

        // Simulation control parameters
        SimType simType_ = SimType::none;
        unsigned int verbosity_ = 0;
        bool strictInput_ = false;
        std::vector<utils::DeckKeyReport> unusedKeys_;
        std::vector<utils::DeckKeyReport> ignoredKeys_;
        unsigned long nTrial_ = 1;
        OutputMode outputMode_ = OutputMode::h5;
        std::string modelName_ = "slug_sim";
        std::string outDir_;
        std::vector<double> outTimes_;
        pdfs::PDF outTimeDist_;
        double intRelTol_ = 1e-2;
        double intAbsTol_ = 1e-3;
        std::size_t intMaxIter_ = 1UL << 19;
        double z_ = 0.0;
        unsigned long checkpointInterval_ = 0;
        std::string inputDeckStr_;
        // Output content flags, parsed by readOutput(); see
        // writeCluster()'s own comment and OutputManager's own
        // constructor for where these are enforced/consumed.
        bool writeCluster_ = true;
        bool writeClusterSpec_ = true;
        bool writeClusterPhot_ = true;
        bool writeGalaxy_ = true;
        bool writeGalaxySpec_ = true;
        bool writeGalaxyPhot_ = true;
        bool writeClusterYields_ = true;
        bool writeGalaxyYields_ = true;
        // Physics settings
        pdfs::PDF imf_;
        pdfs::PDF cmf_;
        pdfs::PDF fehDist_;
        pdfs::PDF sfr_;
        pdfs::PDF sfrDist_;
        pdfs::PDF clf_;
        double fCluster_ = 1.0;
        std::shared_ptr<tracks::Tracks3D> tracks_;
        std::shared_ptr<tracks::Tracks2D> constFeHTracks_;
        double minStochMass_ = 0.0;
        double fracStochMass_ = 1.0;
        std::shared_ptr<specsyn::Specsyn> specsyn_;
        std::shared_ptr<phot::FilterCollection> filters_;
        bool computeLbol_ = false;
        pdfs::PDF avDist_;
        pdfs::PDF avDistField_;
        pdfs::PDF avNebFac_;
        std::shared_ptr<extinct::Extinct> extinct_;
        nebular::NebularControls nebControls_;
        std::shared_ptr<nebular::Nebular> nebular_;
        std::vector<yields::YieldChannelDescriptor> yieldChannels_;
        std::shared_ptr<yields::Yields> yields_;
        bool yieldsChannelDecomposed_ = true;
        bool noDecay_ = false;
        feedback::WRwindModel wrWindModel_ = feedback::WRwindModel::nugisLamers00_;
        feedback::OBwindModel obWindModel_ = feedback::OBwindModel::vinkSander21_;
        feedback::AGBwindModel agbWindModel_ = feedback::AGBwindModel::slug2_;
        feedback::OtherwindModel otherWindModel_ = feedback::OtherwindModel::vesc_;
        // writeClusterFeedback_/writeGalaxyFeedback_ are output content
        // flags like writeCluster_ etc. above, parsed by readOutput(),
        // but declared here instead: the eight above exactly fill one
        // 8-byte word, while these two exactly fill the 2-byte gap left
        // by the six 1-byte members just above before the next double,
        // so declaring them alongside the others would add 8 bytes of
        // padding (and trip clang-analyzer-optin.performance.Padding)
        bool writeClusterFeedback_ = true;
        bool writeGalaxyFeedback_ = true;
        double minIsotopeLifetime_ = yields::defaultMinIsotopeLifetime;
        std::vector<double> snMassLimits_;
        std::shared_ptr<feedback::Winds> winds_;
        // Output wavelength grid (spectra.wl_min, spectra.wl_max,
        // spectra.nwl), read by readSpectra and passed through to
        // whichever spectral synthesizer constructor it builds. Any
        // of the three not supplied in the deck falls back to
        // specsyn::defaultWlMin/defaultWlMax/defaultNWl (see
        // readSpectra's own comment for why); these in-class
        // initializers are only ever visible if readSpectra returns
        // early because spectra.model itself was absent, in which
        // case no spectral synthesizer exists and these values are
        // never read.
        double wlMin_ = 0.0;
        double wlMax_ = 0.0;
        unsigned long nWl_ = 0;
    };

} // namespace io

#endif // SIMCONTROLS_HPP