engine_multiscale.h

Go to the documentation of this file.

#pragma once
 
#include "gridfire/engine/engine_abstract.h"
#include "gridfire/engine/views/engine_view_abstract.h"
#include "gridfire/engine/engine_graph.h"
 
#include "unsupported/Eigen/NonLinearOptimization"
 
namespace gridfire {
    struct QSECacheConfig {
        double T9_tol; 
        double rho_tol; 
        double Yi_tol; 
    };

    struct QSECacheKey {
        double m_T9;
        double m_rho;
        std::vector<double> m_Y; 
 
        std::size_t m_hash = 0; 
 
        // TODO: We should probably sort out how to adjust these from absolute to relative tolerances.
        QSECacheConfig m_cacheConfig = {
            1e-3, // Default tolerance for T9
            1e-1, // Default tolerance for rho
            1e-3  // Default tolerance for species abundances
        };

        QSECacheKey(
            const double T9,
            const double rho,
            const std::vector<double>& Y
        );

        size_t hash() const;

        static long bin(double value, double tol);

        bool operator==(const QSECacheKey& other) const;
 
    };
}
 
// Needs to be in this order (splitting gridfire namespace up) to avoid some issues with forward declarations and the () operator.
namespace std {
    template <>
    struct hash<gridfire::QSECacheKey> {
        size_t operator()(const gridfire::QSECacheKey& key) const noexcept {
            // The hash is pre-computed, so we just return it.
            return key.m_hash;
        }
    };
} // namespace std
 
namespace gridfire {
    class MultiscalePartitioningEngineView final: public DynamicEngine, public EngineView<DynamicEngine> {
        typedef std::tuple<std::vector<fourdst::atomic::Species>, std::vector<size_t>, std::vector<fourdst::atomic::Species>, std::vector<size_t>> QSEPartition;
    public:
        explicit MultiscalePartitioningEngineView(GraphEngine& baseEngine);

        [[nodiscard]] const std::vector<fourdst::atomic::Species> & getNetworkSpecies() const override;

        [[nodiscard]] std::expected<StepDerivatives<double>, expectations::StaleEngineError> calculateRHSAndEnergy(
            const std::vector<double> &Y_full,
            double T9,
            double rho
        ) const override;

        void generateJacobianMatrix(
            const std::vector<double> &Y_full,
            double T9,
            double rho
        ) const override;

        [[nodiscard]] double getJacobianMatrixEntry(
            int i_full,
            int j_full
        ) const override;

        void generateStoichiometryMatrix() override;

        [[nodiscard]] int getStoichiometryMatrixEntry(
            int speciesIndex,
            int reactionIndex
        ) const override;

        [[nodiscard]] double calculateMolarReactionFlow(
            const reaction::Reaction &reaction,
            const std::vector<double> &Y_full,
            double T9,
            double rho
        ) const override;

        [[nodiscard]] const reaction::LogicalReactionSet & getNetworkReactions() const override;

        void setNetworkReactions(
            const reaction::LogicalReactionSet &reactions
        ) override;

        [[nodiscard]] std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> getSpeciesTimescales(
            const std::vector<double> &Y,
            double T9,
            double rho
        ) const override;

        [[nodiscard]] std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> getSpeciesDestructionTimescales(
            const std::vector<double> &Y,
            double T9,
            double rho
        ) const override;

        fourdst::composition::Composition update(
            const NetIn &netIn
        ) override;

        bool isStale(const NetIn& netIn) override;

        void setScreeningModel(
            screening::ScreeningType model
        ) override;

        [[nodiscard]] screening::ScreeningType getScreeningModel() const override;

        const DynamicEngine & getBaseEngine() const override;

        std::vector<std::vector<size_t>> analyzeTimescalePoolConnectivity(
            const std::vector<std::vector<size_t>> &timescale_pools,
            const std::vector<double> &Y,
            double T9,
            double rho
        ) const;

        void partitionNetwork(
            const std::vector<double>& Y,
            double T9,
            double rho
        );

        void partitionNetwork(
            const NetIn& netIn
        );

        void exportToDot(
            const std::string& filename,
            const std::vector<double>& Y,
            const double T9,
            const double rho
        ) const;

        [[nodiscard]] int getSpeciesIndex(const fourdst::atomic::Species &species) const override;

        [[nodiscard]] std::vector<double> mapNetInToMolarAbundanceVector(const NetIn &netIn) const override;

        [[nodiscard]] PrimingReport primeEngine(const NetIn &netIn) override;

        [[nodiscard]] std::vector<fourdst::atomic::Species> getFastSpecies() const;
        [[nodiscard]] const std::vector<fourdst::atomic::Species>& getDynamicSpecies() const;

        fourdst::composition::Composition equilibrateNetwork(
            const std::vector<double> &Y,
            double T9,
            double rho
        );

        fourdst::composition::Composition equilibrateNetwork(
            const NetIn &netIn
        );
 
 
    private:
        struct QSEGroup {
            std::set<size_t> species_indices; 
            bool is_in_equilibrium = false;      
            std::set<size_t> algebraic_indices; 
            std::set<size_t> seed_indices; 
            double mean_timescale; 

            bool operator<(const QSEGroup& other) const;
            bool operator>(const QSEGroup& other) const;
            bool operator==(const QSEGroup& other) const;
            bool operator!=(const QSEGroup& other) const;
        };

        struct EigenFunctor {
            using InputType = Eigen::Matrix<double, Eigen::Dynamic, 1>;
            using OutputType = Eigen::Matrix<double, Eigen::Dynamic, 1>;
            using JacobianType = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>;
            enum {
                InputsAtCompileTime = Eigen::Dynamic,
                ValuesAtCompileTime = Eigen::Dynamic
            };

            MultiscalePartitioningEngineView* m_view;
            const std::vector<size_t>& m_qse_solve_indices;
            const std::vector<double>& m_Y_full_initial;
            const double m_T9;
            const double m_rho;
            const Eigen::VectorXd& m_Y_scale;

            EigenFunctor(
                MultiscalePartitioningEngineView& view,
                const std::vector<size_t>& qse_solve_indices,
                const std::vector<double>& Y_full_initial,
                const double T9,
                const double rho,
                const Eigen::VectorXd& Y_scale
            ) :
            m_view(&view),
            m_qse_solve_indices(qse_solve_indices),
            m_Y_full_initial(Y_full_initial),
            m_T9(T9),
            m_rho(rho),
            m_Y_scale(Y_scale) {}

            [[nodiscard]] int values() const { return m_qse_solve_indices.size(); }
            [[nodiscard]] int inputs() const { return m_qse_solve_indices.size(); }

            int operator()(const InputType& v_qse, OutputType& f_qse) const;
            int df(const InputType& v_qse, JacobianType& J_qse) const;
        };
 

        struct CacheStats {
            enum class operators {
                CalculateRHSAndEnergy,
                GenerateJacobianMatrix,
                CalculateMolarReactionFlow,
                GetSpeciesTimescales,
                GetSpeciesDestructionTimescales,
                Other,
                All
            };

            std::map<operators, std::string> operatorsNameMap = {
                {operators::CalculateRHSAndEnergy, "calculateRHSAndEnergy"},
                {operators::GenerateJacobianMatrix, "generateJacobianMatrix"},
                {operators::CalculateMolarReactionFlow, "calculateMolarReactionFlow"},
                {operators::GetSpeciesTimescales, "getSpeciesTimescales"},
                {operators::GetSpeciesDestructionTimescales, "getSpeciesDestructionTimescales"},
                {operators::Other, "other"}
            };

            size_t m_hit = 0;
            size_t m_miss = 0;

            std::map<operators, size_t> m_operatorHits = {
                {operators::CalculateRHSAndEnergy, 0},
                {operators::GenerateJacobianMatrix, 0},
                {operators::CalculateMolarReactionFlow, 0},
                {operators::GetSpeciesTimescales, 0},
                {operators::GetSpeciesDestructionTimescales, 0},
                {operators::Other, 0}
            };

            std::map<operators, size_t> m_operatorMisses = {
                {operators::CalculateRHSAndEnergy, 0},
                {operators::GenerateJacobianMatrix, 0},
                {operators::CalculateMolarReactionFlow, 0},
                {operators::GetSpeciesTimescales, 0},
                {operators::GetSpeciesDestructionTimescales, 0},
                {operators::Other, 0}
            };

            void hit(const operators op=operators::Other);
            void miss(const operators op=operators::Other);

            [[nodiscard]] size_t hits(const operators op=operators::All) const;
            [[nodiscard]] size_t misses(const operators op=operators::All) const;
        };
 
 
    private:
        quill::Logger* m_logger = LogManager::getInstance().getLogger("log");
        GraphEngine& m_baseEngine;
        std::vector<QSEGroup> m_qse_groups;
        std::vector<fourdst::atomic::Species> m_dynamic_species;
        std::vector<size_t> m_dynamic_species_indices;
        std::vector<fourdst::atomic::Species> m_algebraic_species;
        std::vector<size_t> m_algebraic_species_indices;

        std::vector<size_t> m_activeSpeciesIndices;
        std::vector<size_t> m_activeReactionIndices;
 
        // TODO: Enhance the hashing for the cache to consider not just T and rho but also the current abundance in some careful way that automatically ignores small changes (i.e. network should only be repartitioned sometimes)
        mutable std::unordered_map<QSECacheKey, std::vector<double>> m_qse_abundance_cache;
        mutable CacheStats m_cacheStats;
 
 
 
 
    private:
        std::vector<std::vector<size_t>> partitionByTimescale(
            const std::vector<double> &Y_full,
            double T9,
            double rho
        ) const;

        std::unordered_map<size_t, std::vector<size_t>> buildConnectivityGraph(
            const std::unordered_set<size_t> &fast_reaction_indices
        ) const;

        std::vector<QSEGroup> validateGroupsWithFluxAnalysis(
            const std::vector<QSEGroup> &candidate_groups,
            const std::vector<double>& Y,
            double T9,
            double rho
        ) const;

        std::vector<double> solveQSEAbundances(
            const std::vector<double> &Y_full,
            double T9,
            double rho
        );

        size_t identifyMeanSlowestPool(
            const std::vector<std::vector<size_t>>& pools,
            const std::vector<double> &Y,
            double T9,
            double rho
        ) const;

        std::unordered_map<size_t, std::vector<size_t>> buildConnectivityGraph(
            const std::vector<size_t>& species_pool
        ) const;

        std::vector<QSEGroup> constructCandidateGroups(
            const std::vector<std::vector<size_t>>& candidate_pools,
            const std::vector<double>& Y,
            double T9,
            double rho
        ) const;
    };
}