perf(graph_engine): finished sparsity system for jacobian, major preformance win, roughly 20x faster

essentially all callers can now inform the graph engine about which species they hold active and graph engine then uses those to define a sparsity pattern and only calculate the jacobian along that sparsity pattern

src/include/gridfire/engine/engine_abstract.h

@@ -157,14 +157,19 @@ namespace gridfire {
             double rho
         ) const = 0;
 
+        virtual void generateJacobianMatrix(
+            const fourdst::composition::Composition& comp,
+            double T9,
+            double rho,
+            const std::vector<fourdst::atomic::Species>& activeSpecies
+        ) const = 0;
+
         virtual void generateJacobianMatrix(
             const fourdst::composition::Composition& comp,
             double T9,
             double rho,
             const SparsityPattern& sparsityPattern
-        ) const {
-            throw std::logic_error("Sparsity pattern not supported by this engine.");
-        }
+        ) const = 0;
 
         /**
          * @brief Get an entry from the previously generated Jacobian matrix.

src/include/gridfire/engine/engine_graph.h

@@ -12,7 +12,6 @@
 #include "gridfire/screening/screening_types.h"
 #include "gridfire/partition/partition_abstract.h"
 #include "gridfire/engine/procedures/construction.h"
-#include "gridfire/utils/general_composition.h"
 
 #include <string>
 #include <unordered_map>
@@ -235,6 +234,25 @@ namespace gridfire {
             double rho
         ) const override;
 
+        /**
+         * @brief Generates the Jacobian matrix for the current state with a specified set of active species.
+         *        generally this will be much faster than the full matrix generation. Here we use forward mode
+         *        to generate the Jacobian only for the active species.
+         * @param comp The Composition object containing current abundances.
+         * @param T9 The temperature in units of 10^9 K.
+         * @param rho The density in g/cm^3.
+         * @param activeSpecies A vector of Species objects representing the active species.
+         *
+         * @see getJacobianMatrixEntry()
+         * @see generateJacobianMatrix()
+         */
+        void generateJacobianMatrix(
+            const fourdst::composition::Composition& comp,
+            double T9,
+            double rho,
+            const std::vector<fourdst::atomic::Species>& activeSpecies
+        ) const override;
+
         /**
          * @brief Generates the Jacobian matrix for the current state with a specified sparsity pattern.
          *
@@ -717,6 +735,7 @@ namespace gridfire {
             // Forward cacheing
             size_t reaction_index;
             reaction::ReactionType reaction_type;
+            uint64_t reaction_hash;
             std::vector<size_t> unique_reactant_indices;
             std::vector<int> reactant_powers;
             double symmetry_factor;
@@ -800,6 +819,7 @@ namespace gridfire {
         mutable CppAD::ADFun<double> m_epsADFun; ///< CppAD function for the energy generation rate.
         mutable CppAD::sparse_jac_work m_jac_work; ///< Work object for sparse Jacobian calculations.
         CppAD::sparse_rc<std::vector<size_t>> m_full_jacobian_sparsity_pattern; ///< Full sparsity pattern for the Jacobian matrix.
+        std::set<std::pair<size_t, size_t>> m_full_sparsity_set; ///< For quick lookups of the base sparsity pattern
 
         std::vector<std::unique_ptr<AtomicReverseRate>> m_atomicReverseRates;
 
@@ -813,6 +833,7 @@ namespace gridfire {
         BuildDepthType m_depth;
 
         std::vector<PrecomputedReaction> m_precomputedReactions; ///< Precomputed reactions for efficiency.
+        std::unordered_map<uint64_t, size_t> m_precomputedReactionIndexMap; ///< Set of hashed precomputed reactions for quick lookup.
         std::unique_ptr<partition::PartitionFunction> m_partitionFunction; ///< Partition function for the network.
 
     private:
@@ -890,10 +911,10 @@ namespace gridfire {
 
         [[nodiscard]] StepDerivatives<double> calculateAllDerivativesUsingPrecomputation(
             const fourdst::composition::Composition &comp,
-            const std::vector<double>& bare_rates,
+            const std::vector<double> &bare_rates,
             const std::vector<double> &bare_reverse_rates,
             double T9,
-            double rho
+            double rho, const reaction::ReactionSet &activeReactions
         ) const;
 
         /**

src/include/gridfire/engine/views/engine_adaptive.h

@@ -140,6 +140,20 @@ namespace gridfire {
             double rho
         ) const override;
 
+        void generateJacobianMatrix(
+            const fourdst::composition::Composition &comp,
+            double T9,
+            double rho,
+            const std::vector<fourdst::atomic::Species> &activeSpecies
+        ) const override;
+
+        void generateJacobianMatrix(
+            const fourdst::composition::Composition &comp,
+            double T9,
+            double rho,
+            const SparsityPattern &sparsityPattern
+        ) const override;
+
         /**
          * @brief Gets an entry from the Jacobian matrix for the active species.
          *
@@ -409,7 +423,7 @@ namespace gridfire {
          * 5. For each reaction, it calls the base engine's `calculateMolarReactionFlow` to get the flow rate.
          * 6. Stores the reaction pointer and its flow rate in a `ReactionFlow` struct and adds it to the returned vector.
          */
-        std::pair<std::vector<ReactionFlow>, fourdst::composition::Composition> calculateAllReactionFlows(
+        [[nodiscard]] std::pair<std::vector<ReactionFlow>, fourdst::composition::Composition> calculateAllReactionFlows(
             const NetIn& netIn
         ) const;
         /**

src/include/gridfire/engine/views/engine_defined.h

@@ -65,9 +65,44 @@ namespace gridfire{
          */
         void generateJacobianMatrix(
             const fourdst::composition::Composition& comp,
-            const double T9,
-            const double rho
+            double T9,
+            double rho
         ) const override;
+
+        /**
+         * @brief Generates the Jacobian matrix for the active species.
+         *
+         * @param comp A Composition object containing the current composition of the system
+         * @param T9 The temperature in units of 10^9 K.
+         * @param rho The density in g/cm^3.
+         * @param activeSpecies The vector of active species to include in the Jacobian.
+         *
+         * @throws std::runtime_error If the view is stale.
+         */
+        void generateJacobianMatrix(
+            const fourdst::composition::Composition &comp,
+            double T9,
+            double rho,
+            const std::vector<fourdst::atomic::Species> &activeSpecies
+        ) const override;
+
+        /**
+         * @brief Generates the Jacobian matrix for a given sparsity pattern
+         *
+         * @param comp A Composition object containing the current composition of the system
+         * @param T9 The temperature in units of 10^9 K.
+         * @param rho The density in g/cm^3.
+         * @param sparsityPattern The sparsity pattern to use for the Jacobian matrix.
+         *
+         * @throws std::runtime_error If the view is stale.
+         */
+        void generateJacobianMatrix(
+            const fourdst::composition::Composition &comp,
+            double T9,
+            double rho,
+            const SparsityPattern &sparsityPattern
+        ) const override;
+
         /**
          * @brief Gets an entry from the Jacobian matrix for the active species.
          *

src/include/gridfire/engine/views/engine_multiscale.h

@@ -271,6 +271,64 @@ namespace gridfire {
             double rho
         ) const override;
 
+        /**
+         * @brief Generates the Jacobian matrix for a subset of active species.
+         *
+         * @param comp The current composition.
+         * @param T9 Temperature in units of 10^9 K.
+         * @param rho Density in g/cm^3.
+         * @param activeSpecies The subset of species to include in the Jacobian.
+         *
+         * @par Purpose
+         * To compute a reduced Jacobian matrix for implicit solvers that only
+         * consider a subset of species.
+         *
+         * @par How
+         * Similar to the full Jacobian generation, it first checks the QSE cache.
+         * On a hit, it calls the base engine's `generateJacobianMatrix` with
+         * the specified active species. The returned Jacobian still reflects
+         * the full network, but only for the active species subset.
+         *
+         * @pre The engine must have a valid QSE cache entry for the given state.
+         * @post The base engine's internal Jacobian is updated for the active species.
+         *
+         * @throws exceptions::StaleEngineError If the QSE cache misses.
+         */
+        void generateJacobianMatrix(
+            const fourdst::composition::Composition &comp,
+            double T9,
+            double rho,
+            const std::vector<fourdst::atomic::Species> &activeSpecies
+        ) const override;
+
+        /**
+         * @brief Generates the Jacobian matrix using a sparsity pattern.
+         *
+         * @param comp The current composition.
+         * @param T9 Temperature in units of 10^9 K.
+         * @param rho Density in g/cm^3.
+         * @param sparsityPattern The sparsity pattern to use for the Jacobian.
+         *
+         * @par Purpose
+         * To compute the Jacobian matrix while leveraging a known sparsity pattern
+         * for efficiency. This is effectively a lower level version of the active species method.
+         *
+         * @par How
+         * It first checks the QSE cache. On a hit, it delegates to the base engine's
+         * `generateJacobianMatrix` method with the provided sparsity pattern.
+         *
+         * @pre The engine must have a valid QSE cache entry for the given state.
+         * @post The base engine's internal Jacobian is updated according to the sparsity pattern.
+         *
+         * @throws exceptions::StaleEngineError If the QSE cache misses.
+         */
+        void generateJacobianMatrix(
+            const fourdst::composition::Composition &comp,
+            double T9,
+            double rho,
+            const SparsityPattern &sparsityPattern
+        ) const override;
+
         /**
          * @brief Gets an entry from the previously generated Jacobian matrix.
          *
@@ -720,6 +778,12 @@ namespace gridfire {
             const NetIn &netIn
         );
 
+        bool involvesSpecies(const fourdst::atomic::Species &species) const;
+
+        bool involvesSpeciesInQSE(const fourdst::atomic::Species &species) const;
+
+        bool involvesSpeciesInDynamic(const fourdst::atomic::Species &species) const;
+
 
     private:
         /**
@@ -806,7 +870,7 @@ namespace gridfire {
              */
             MultiscalePartitioningEngineView* m_view;
             /**
-             * @brief Indices of the species to solve for in the QSE group.
+             * @brief The set of species to solve for in the QSE group.
              */
             const std::set<fourdst::atomic::Species>& m_qse_solve_species;
             /**

src/include/gridfire/engine/views/engine_priming.h

@@ -1,9 +1,7 @@
 #pragma once
 
-#include "gridfire/engine/engine_abstract.h"
 #include "gridfire/engine/views/engine_defined.h"
 
-
 #include "fourdst/logging/logging.h"
 #include "fourdst/composition/atomicSpecies.h"
 
@@ -51,7 +49,7 @@ namespace gridfire {
 
 
     private:
-        quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log");
+        quill::Logger* m_logger = LogManager::getInstance().getLogger("log");
         fourdst::atomic::Species m_primingSpecies; ///< The priming species, if specified.
     private:
         /**

src/include/gridfire/exceptions/error_utils.h

@@ -0,0 +1,23 @@
+#pragma once
+
+#include <exception>
+#include <string>
+#include <utility>
+
+namespace gridfire::exceptions {
+    class UtilityError : public std::exception {};
+
+    class HashingError final : public UtilityError {
+    public:
+        explicit HashingError() = default;
+
+        explicit HashingError(std::string message)
+            : m_message(std::move(message)) {}
+
+        [[nodiscard]] const char* what() const noexcept override {
+            return m_message.c_str();
+        }
+    private:
+        std::string m_message;
+    };
+}

src/include/gridfire/exceptions/exceptions.h

@@ -1,3 +1,4 @@
 #pragma once
 
 #include "gridfire/exceptions/error_engine.h"
+#include "gridfire/exceptions/error_utils.h"

src/include/gridfire/utils/hashing.h

@@ -2,6 +2,9 @@
 
 #include <cstdint>
 
+#include "gridfire/exceptions/exceptions.h"
+#include "gridfire/reaction/reaction.h"
+
 namespace gridfire::utils {
     /**
      * @brief Generate a unique hash for an isotope given its mass number (A) and atomic number (Z).
@@ -16,4 +19,47 @@ namespace gridfire::utils {
         return (static_cast<uint_fast32_t>(a) << 8) | static_cast<uint_fast32_t>(z);
     }
 
+    namespace hashing::reaction {
+        static std::uint64_t splitmix64(std::uint64_t x) noexcept {
+            x += 0x9E3779B97F4A7C15ULL;
+            x = (x ^ (x >> 30)) * 0xBF58476D1CE4E5B9ULL;
+            x = (x ^ (x >> 27)) * 0x94D049BB133111EBULL;
+            x ^= (x >> 31);
+            return x;
+        }
+
+        static std::uint64_t mix_species(const unsigned a, const unsigned z) noexcept {
+            const std::uint64_t code = (static_cast<std::uint64_t>(a) << 7) | static_cast<std::uint64_t>(z);
+            return splitmix64(code);
+        }
+
+        static std::uint64_t multiset_combine(std::uint64_t acc, const std::uint64_t x) noexcept {
+            acc += x;
+            acc ^= (x << 23) | (x >> (64 - 23));
+            acc = splitmix64(acc);
+            return acc;
+        }
+    }
+
+    inline std::uint64_t hash_reaction(const reaction::Reaction& reaction) noexcept {
+        using namespace hashing::reaction;
+
+        std::uint64_t hR = 0;
+
+        for (const auto& s : reaction.reactants()) {
+            hR = multiset_combine(hR, mix_species(static_cast<unsigned>(s.a()),
+                                                  static_cast<unsigned>(s.z())));
+        }
+
+        std::uint64_t hP = 0;
+        for (const auto& s : reaction.products()) {
+            hP = multiset_combine(hP, mix_species(static_cast<unsigned>(s.a()),
+                                                  static_cast<unsigned>(s.z())));
+        }
+
+        std::uint64_t h = splitmix64(hR ^ 0xC3A5C85C97CB3127ULL);
+        h ^= splitmix64((hP << 1) | (hP >> 63));
+        return splitmix64(h);
+    }
+
 }

src/lib/engine/engine_graph.cpp

@@ -5,6 +5,7 @@
 #include "gridfire/engine/procedures/priming.h"
 #include "gridfire/partition/partition_ground.h"
 #include "gridfire/engine/procedures/construction.h"
+#include "gridfire/utils/hashing.h"
 
 #include "fourdst/composition/species.h"
 #include "fourdst/composition/atomicSpecies.h"
@@ -75,12 +76,11 @@ namespace gridfire {
         if (m_usePrecomputation) {
             std::vector<double> bare_rates;
             std::vector<double> bare_reverse_rates;
-            bare_rates.reserve(m_reactions.size());
-            bare_reverse_rates.reserve(m_reactions.size());
+            bare_rates.reserve(activeReactions.size());
+            bare_reverse_rates.reserve(activeReactions.size());
 
-            // TODO: Add cache to this
-
-            for (const auto& reaction: m_reactions) {
+            for (const auto& reaction: activeReactions) {
+                assert(m_reactions.contains(*reaction)); // A bug which results in this failing indicates a serious internal inconsistency and should only be present during development.
                 bare_rates.push_back(reaction->calculate_rate(T9, rho, Ye, mue, comp.getMolarAbundanceVector(), m_indexToSpeciesMap));
                 if (reaction->type() != reaction::ReactionType::WEAK) {
                     bare_reverse_rates.push_back(calculateReverseRate(*reaction, T9, rho, comp));
@@ -88,7 +88,7 @@ namespace gridfire {
             }
 
             // --- The public facing interface can always use the precomputed version since taping is done internally ---
-            return calculateAllDerivativesUsingPrecomputation(comp, bare_rates, bare_reverse_rates, T9, rho);
+            return calculateAllDerivativesUsingPrecomputation(comp, bare_rates, bare_reverse_rates, T9, rho, activeReactions);
         } else {
             return calculateAllDerivatives<double>(
                 comp.getMolarAbundanceVector(),
@@ -171,6 +171,8 @@ namespace gridfire {
         collectAtomicReverseRateAtomicBases();
         generateStoichiometryMatrix();
         reserveJacobianMatrix();
+
+        // PERF: These do *a lot* of the same work* can this be optimized to only do the common work once?
         recordADTape(); // Record the AD tape for the RHS function
         recordEpsADTape(); // Record the AD tape for the energy generation rate function
 
@@ -183,6 +185,17 @@ namespace gridfire {
         m_rhsADFun.subgraph_sparsity(select_domain, select_range, false, m_full_jacobian_sparsity_pattern);
         m_jac_work.clear();
 
+        m_full_sparsity_set.clear();
+        const auto& rows = m_full_jacobian_sparsity_pattern.row();
+        const auto& cols = m_full_jacobian_sparsity_pattern.col();
+        const size_t nnz = m_full_jacobian_sparsity_pattern.nnz();
+
+        for (size_t k = 0; k < nnz; ++k) {
+            if (cols[k] < m_networkSpecies.size()) {
+                m_full_sparsity_set.insert(std::make_pair(rows[k], cols[k]));
+            }
+        }
+
         precomputeNetwork();
         LOG_INFO(m_logger, "Internal maps synchronized. Network contains {} species and {} reactions.",
                  m_networkSpecies.size(), m_reactions.size());
@@ -190,7 +203,6 @@ namespace gridfire {
 
     // --- Network Graph Construction Methods ---
     void GraphEngine::collectNetworkSpecies() {
-        //TODO: Ensure consistent ordering in the m_networkSpecies vector so that it is ordered by species mass.
         m_networkSpecies.clear();
         m_networkSpeciesMap.clear();
 
@@ -217,7 +229,7 @@ namespace gridfire {
             }
         }
         // TODO: Currently this works. We sort the vector based on mass so that for the same set of species we always get the same ordering and we get the same ordering as a composition with the same set of species
-        //       However, we need some checks so that when we get a composition we confirm that it is the same ordering / contains teh same species. This is important for the ODE integrator to work properly.
+        //       However, we need some checks so that when we get a composition we confirm that it is the same ordering / contains the same species. This is important for the ODE integrator to work properly.
         std::ranges::sort(m_networkSpecies, [](const fourdst::atomic::Species& a, const fourdst::atomic::Species& b) -> bool {
             return a.mass() < b.mass(); // Otherwise, sort by mass
         });
@@ -255,14 +267,10 @@ namespace gridfire {
 
     // --- Basic Accessors and Queries ---
     const std::vector<fourdst::atomic::Species>& GraphEngine::getNetworkSpecies() const {
-        // Returns a constant reference to the vector of unique species in the network.
-        // LOG_TRACE_L3(m_logger, "Providing access to network species vector. Size: {}.", m_networkSpecies.size());
         return m_networkSpecies;
     }
 
     const reaction::ReactionSet& GraphEngine::getNetworkReactions() const {
-        // Returns a constant reference to the set of reactions in the network.
-        // LOG_TRACE_L3(m_logger, "Providing access to network reactions set. Size: {}.", m_reactions.size());
         return m_reactions;
     }
 
@@ -573,64 +581,75 @@ namespace gridfire {
         const std::vector<double> &bare_rates,
         const std::vector<double> &bare_reverse_rates,
         const double T9,
-        const double rho
+        const double rho,
+        const reaction::ReactionSet &activeReactions
     ) const {
         // --- Calculate screening factors ---
         const std::vector<double> screeningFactors = m_screeningModel->calculateScreeningFactors(
-            m_reactions,
+            activeReactions,
             m_networkSpecies,
             comp.getMolarAbundanceVector(),
             T9,
             rho
         );
 
-        // TODO: Fix up the precomputation to use the new comp in interface as opposed to a raw vector of molar abundances
-        //       This will require carefully checking the way the precomputation is stashed.
-
         // --- Optimized loop ---
         std::vector<double> molarReactionFlows;
         molarReactionFlows.reserve(m_precomputedReactions.size());
 
-        for (const auto& precomp : m_precomputedReactions) {
+        size_t reactionCounter = 0;
+        for (const auto& reaction : activeReactions) {
+            // --- Efficient lookup of only the active reactions ---
+            uint64_t reactionHash = utils::hash_reaction(*reaction);
+            const size_t reactionIndex = m_precomputedReactionIndexMap.at(reactionHash);
+            PrecomputedReaction precomputedReaction = m_precomputedReactions[reactionIndex];
+
+            // --- Forward abundance product ---
             double forwardAbundanceProduct = 1.0;
-            for (size_t i = 0; i < precomp.unique_reactant_indices.size(); ++i) {
-                const size_t reactantIndex = precomp.unique_reactant_indices[i];
+            for (size_t i = 0; i < precomputedReaction.unique_reactant_indices.size(); ++i) {
+                const size_t reactantIndex = precomputedReaction.unique_reactant_indices[i];
                 const fourdst::atomic::Species& reactant = m_networkSpecies[reactantIndex];
-                const int power = precomp.reactant_powers[i];
+                const int power = precomputedReaction.reactant_powers[i];
 
                 forwardAbundanceProduct *= std::pow(comp.getMolarAbundance(reactant), power);
             }
 
-            const double bare_rate = bare_rates[precomp.reaction_index];
-            const double screeningFactor = screeningFactors[precomp.reaction_index];
-            const size_t numReactants = m_reactions[precomp.reaction_index].reactants().size();
-            const size_t numProducts = m_reactions[precomp.reaction_index].products().size();
+            const double bare_rate = bare_rates.at(reactionCounter);
 
+            const double screeningFactor = screeningFactors[reactionCounter];
+            const size_t numReactants = m_reactions[reactionIndex].reactants().size();
+            const size_t numProducts = m_reactions[reactionIndex].products().size();
+
+            // --- Forward reaction flow ---
             const double forwardMolarReactionFlow =
                     screeningFactor *
                     bare_rate *
-                    precomp.symmetry_factor *
+                    precomputedReaction.symmetry_factor *
                     forwardAbundanceProduct *
                     std::pow(rho, numReactants >  1 ? static_cast<double>(numReactants) - 1 : 0.0);
 
+            // --- Reverse reaction flow ---
+            // Only do this is the reaction has a non-zero reverse symmetry factor (i.e. is reversible)
             double reverseMolarReactionFlow = 0.0;
-            if (precomp.reverse_symmetry_factor != 0.0 and m_useReverseReactions) {
-                const double bare_reverse_rate = bare_reverse_rates[precomp.reaction_index];
+            if (precomputedReaction.reverse_symmetry_factor != 0.0 and m_useReverseReactions) {
+                const double bare_reverse_rate = bare_reverse_rates.at(reactionCounter);
+
                 double reverseAbundanceProduct = 1.0;
-                for (size_t i = 0; i < precomp.unique_product_indices.size(); ++i) {
-                    const size_t productIndex = precomp.unique_product_indices[i];
+                for (size_t i = 0; i < precomputedReaction.unique_product_indices.size(); ++i) {
+                    const size_t productIndex = precomputedReaction.unique_product_indices[i];
                     const fourdst::atomic::Species& product = m_networkSpecies[productIndex];
-                    reverseAbundanceProduct *= std::pow(comp.getMolarAbundance(product), precomp.product_powers[i]);
+                    reverseAbundanceProduct *= std::pow(comp.getMolarAbundance(product), precomputedReaction.product_powers[i]);
                 }
+
                 reverseMolarReactionFlow = screeningFactor *
                     bare_reverse_rate *
-                    precomp.reverse_symmetry_factor *
+                    precomputedReaction.reverse_symmetry_factor *
                     reverseAbundanceProduct *
                     std::pow(rho, numProducts > 1 ? static_cast<double>(numProducts) - 1 : 0.0);
             }
 
             molarReactionFlows.push_back(forwardMolarReactionFlow - reverseMolarReactionFlow);
-
+            reactionCounter++;
         }
 
         // --- Assemble molar abundance derivatives ---
@@ -638,9 +657,12 @@ namespace gridfire {
         for (const auto& species: m_networkSpecies) {
             result.dydt[species] = 0.0; // Initialize the change in abundance for each network species to 0
         }
-        for (size_t j = 0; j < m_precomputedReactions.size(); ++j) {
+
+        reactionCounter = 0;
+        for (const auto& reaction: activeReactions) {
+            size_t j = m_precomputedReactionIndexMap.at(utils::hash_reaction(*reaction));
             const auto& precomp = m_precomputedReactions[j];
-            const double R_j = molarReactionFlows[j];
+            const double R_j = molarReactionFlows[reactionCounter];
 
             for (size_t i = 0; i < precomp.affected_species_indices.size(); ++i) {
                 const size_t speciesIndex = precomp.affected_species_indices[i];
@@ -651,6 +673,7 @@ namespace gridfire {
                 // Update the derivative for this species
                 result.dydt.at(species) += static_cast<double>(stoichiometricCoefficient) * R_j;
             }
+            reactionCounter++;
         }
 
         // --- Calculate the nuclear energy generation rate ---
@@ -802,12 +825,51 @@ namespace gridfire {
         LOG_TRACE_L1_LIMIT_EVERY_N(1000, m_logger, "Jacobian matrix generated with dimensions: {} rows x {} columns.", m_jacobianMatrix.size1(), m_jacobianMatrix.size2());
     }
 
+    void GraphEngine::generateJacobianMatrix(
+        const fourdst::composition::Composition &comp,
+        const double T9,
+        const double rho,
+        const std::vector<fourdst::atomic::Species> &activeSpecies
+    ) const {
+        // PERF: For small k it may make sense to implement a purley forward mode AD computation, some heuristic could be used to switch between the two methods based on k and total network species
+        const size_t k_active = activeSpecies.size();
+
+        // --- 1. Get the list of global indices ---
+        std::vector<size_t> active_indices;
+        active_indices.reserve(k_active);
+
+        for (const auto& species : activeSpecies) {
+            assert(involvesSpecies(species));
+            active_indices.push_back(getSpeciesIndex(species));
+        }
+
+        // --- 2. Build the k x k sparsity pattern ---
+        SparsityPattern sparsityPattern;
+        sparsityPattern.reserve(k_active * k_active);
+
+        for (const size_t i_global : active_indices) { // k rows
+            for (const size_t j_global : active_indices) { // k columns
+                sparsityPattern.emplace_back(i_global, j_global);
+            }
+        }
+
+        // --- 3. Call the sparse reverse-mode implementation ---
+        generateJacobianMatrix(comp, T9, rho, sparsityPattern);
+    }
+
     void GraphEngine::generateJacobianMatrix(
         const fourdst::composition::Composition &comp,
         const double T9,
         const double rho,
         const SparsityPattern &sparsityPattern
     ) const {
+        SparsityPattern intersectionSparsityPattern;
+        for (const auto& entry : sparsityPattern) {
+            if (m_full_sparsity_set.contains(entry)) {
+                intersectionSparsityPattern.push_back(entry);
+            }
+        }
+
         // --- Pack the input variables into a vector for CppAD ---
         const size_t numSpecies = m_networkSpecies.size();
         std::vector<double> x(numSpecies + 2, 0.0);
@@ -819,13 +881,13 @@ namespace gridfire {
         x[numSpecies + 1] = rho;
 
         // --- Convert into CppAD Sparsity pattern ---
-        const size_t nnz = sparsityPattern.size(); // Number of non-zero entries in the sparsity pattern
+        const size_t nnz = intersectionSparsityPattern.size(); // Number of non-zero entries in the sparsity pattern
         std::vector<size_t> row_indices(nnz);
         std::vector<size_t> col_indices(nnz);
 
         for (size_t k = 0; k < nnz; ++k) {
-            row_indices[k] = sparsityPattern[k].first;
-            col_indices[k] = sparsityPattern[k].second;
+            row_indices[k] = intersectionSparsityPattern[k].first;
+            col_indices[k] = intersectionSparsityPattern[k].second;
         }
 
         std::vector<double> values(nnz);
@@ -834,7 +896,7 @@ namespace gridfire {
 
         CppAD::sparse_rc<std::vector<size_t>> CppAD_sparsity_pattern(num_rows_jac, num_cols_jac, nnz);
         for (size_t k = 0; k < nnz; ++k) {
-            CppAD_sparsity_pattern.set(k, sparsityPattern[k].first, sparsityPattern[k].second);
+            CppAD_sparsity_pattern.set(k, intersectionSparsityPattern[k].first, intersectionSparsityPattern[k].second);
         }
 
         CppAD::sparse_rcv<std::vector<size_t>, std::vector<double>> jac_subset(CppAD_sparsity_pattern);
@@ -854,7 +916,7 @@ namespace gridfire {
             const size_t col = jac_subset.col()[k];
             const double value = jac_subset.val()[k];
 
-            if (std::abs(value) > MIN_JACOBIAN_THRESHOLD) {
+            if (std::abs(value) > MIN_JACOBIAN_THRESHOLD || row == col) { // Always keep diagonal elements to avoid pathological stiffness
                 m_jacobianMatrix(row, col) = value; // Insert into the sparse matrix
             }
         }
@@ -1241,12 +1303,18 @@ namespace gridfire {
 
         m_precomputedReactions.clear();
         m_precomputedReactions.reserve(m_reactions.size());
+        m_precomputedReactionIndexMap.clear();
+        m_precomputedReactionIndexMap.reserve(m_reactions.size());
 
         for (size_t i = 0; i < m_reactions.size(); ++i) {
             const auto& reaction = m_reactions[i];
             PrecomputedReaction precomp;
             precomp.reaction_index = i;
             precomp.reaction_type = reaction.type();
+            uint64_t reactionHash = utils::hash_reaction(reaction);
+
+            precomp.reaction_hash = reactionHash;
+            m_precomputedReactionIndexMap[reactionHash] = i;
 
             // --- Precompute forward reaction information ---
             // Count occurrences for each reactant to determine powers and symmetry
@@ -1298,6 +1366,7 @@ namespace gridfire {
 
             m_precomputedReactions.push_back(std::move(precomp));
         }
+        LOG_TRACE_L1(m_logger, "Pre-computation complete. Precomputed data for {} reactions.", m_precomputedReactions.size());
     }
 
     bool GraphEngine::AtomicReverseRate::forward(

src/lib/engine/procedures/priming.cpp

@@ -231,11 +231,11 @@ namespace gridfire {
                     // Store the request instead of applying it immediately.
                     requests.push_back({primingSpecies, equilibriumMassFraction, dominantChannel->reactants()});
                 } else {
-                    LOG_ERROR(logger, "Failed to find dominant creation channel for {}.", primingSpecies.name());
+                    LOG_TRACE_L1(logger, "Failed to find dominant creation channel for {}.", primingSpecies.name());
                     report.status = PrimingReportStatus::FAILED_TO_FIND_CREATION_CHANNEL;
                 }
             } else {
-                LOG_WARNING(logger, "No destruction channel found for {}. Using fallback abundance.", primingSpecies.name());
+                LOG_TRACE_L2(logger, "No destruction channel found for {}. Using fallback abundance.", primingSpecies.name());
                 // For species with no destruction, we can't calculate an equilibrium.
                 // We add a request with a tiny fallback abundance to ensure it exists in the network.
                 requests.push_back({primingSpecies, 1e-40, {}});

src/lib/engine/views/engine_adaptive.cpp

@@ -176,9 +176,29 @@ namespace gridfire {
         const fourdst::composition::Composition &comp,
         const double T9,
         const double rho
+    ) const {
+        generateJacobianMatrix(comp, T9, rho, m_activeSpecies);
+    }
+
+    void AdaptiveEngineView::generateJacobianMatrix(
+        const fourdst::composition::Composition &comp,
+        const double T9,
+        const double rho,
+        const std::vector<Species> &activeSpecies
     ) const {
         validateState();
-        m_baseEngine.generateJacobianMatrix(comp, T9, rho);
+        m_baseEngine.generateJacobianMatrix(comp, T9, rho, activeSpecies);
+
+    }
+
+    void AdaptiveEngineView::generateJacobianMatrix(
+        const fourdst::composition::Composition &comp,
+        const double T9,
+        const double rho,
+        const SparsityPattern &sparsityPattern
+    ) const {
+        validateState();
+        m_baseEngine.generateJacobianMatrix(comp, T9, rho, sparsityPattern);
     }
 
     double AdaptiveEngineView::getJacobianMatrixEntry(

src/lib/engine/views/engine_defined.cpp

@@ -86,11 +86,39 @@ namespace gridfire {
         const double rho
     ) const {
         validateNetworkState();
-
+        if (!m_activeSpeciesVectorCache.has_value()) {
+            m_activeSpeciesVectorCache = std::vector<Species>(m_activeSpecies.begin(), m_activeSpecies.end());
+        }
         const MaskedComposition masked(comp, m_activeSpecies);
+        m_baseEngine.generateJacobianMatrix(masked, T9, rho, m_activeSpeciesVectorCache.value());
+    }
 
-        // TODO: We likely want to be able to think more carefully about this so that the jacobian matches the active species/reactions
-        m_baseEngine.generateJacobianMatrix(masked, T9, rho);
+    void DefinedEngineView::generateJacobianMatrix(
+        const fourdst::composition::Composition &comp,
+        const double T9,
+        const double rho,
+        const std::vector<fourdst::atomic::Species> &activeSpecies
+    ) const {
+        validateNetworkState();
+
+        const std::set<fourdst::atomic::Species> activeSpeciesSet(
+            activeSpecies.begin(),
+            activeSpecies.end()
+        );
+
+        const MaskedComposition masked(comp, activeSpeciesSet);
+        m_baseEngine.generateJacobianMatrix(masked, T9, rho, activeSpecies);
+    }
+
+    void DefinedEngineView::generateJacobianMatrix(
+        const fourdst::composition::Composition &comp,
+        const double T9,
+        const double rho,
+        const SparsityPattern &sparsityPattern
+    ) const {
+        validateNetworkState();
+        const MaskedComposition masked(comp, m_activeSpecies);
+        m_baseEngine.generateJacobianMatrix(masked, T9, rho, sparsityPattern);
     }
 
     double DefinedEngineView::getJacobianMatrixEntry(

src/lib/engine/views/engine_multiscale.cpp

@@ -17,8 +17,6 @@
 #include "quill/LogMacros.h"
 #include "quill/Logger.h"
 
-static std::ofstream debug_multiscale_log("debug_multiscale_log.txt");
-
 namespace {
     using namespace fourdst::atomic;
     //TODO: Replace all calls to this function with composition.getMolarAbundanceVector() so that
@@ -214,8 +212,41 @@ namespace gridfire {
         const double T9,
         const double rho
     ) const {
-        // TODO: Add sparsity pattern to this to prevent base engine from doing unnecessary work.
-        m_baseEngine.generateJacobianMatrix(comp, T9, rho);
+        // We do not need to generate the jacobian for QSE species since those entries are by definition 0
+        m_baseEngine.generateJacobianMatrix(comp, T9, rho, m_dynamic_species);
+    }
+
+    void MultiscalePartitioningEngineView::generateJacobianMatrix(
+        const fourdst::composition::Composition &comp,
+        const double T9,
+        const double rho,
+        const std::vector<Species> &activeSpecies
+    ) const {
+        const bool activeSpeciesIsSubset = std::ranges::any_of(activeSpecies, [&](const auto& species) -> bool {
+            return !involvesSpecies(species);
+        });
+        if (activeSpeciesIsSubset) {
+            LOG_CRITICAL(m_logger, "Active species contains species not in the network partition. Cannot generate Jacobian matrix for active species.");
+            throw std::runtime_error("Active species contains species not in the network partition. Cannot generate Jacobian matrix for active species.");
+        }
+
+        std::vector<Species> dynamicActiveSpeciesIntersection;
+        for (const auto& species : activeSpecies) {
+            if (involvesSpeciesInDynamic(species)) {
+                dynamicActiveSpeciesIntersection.push_back(species);
+            }
+        }
+
+        m_baseEngine.generateJacobianMatrix(comp, T9, rho, dynamicActiveSpeciesIntersection);
+    }
+
+    void MultiscalePartitioningEngineView::generateJacobianMatrix(
+        const fourdst::composition::Composition &comp,
+        const double T9,
+        const double rho,
+        const SparsityPattern &sparsityPattern
+    ) const {
+        return m_baseEngine.generateJacobianMatrix(comp, T9, rho, sparsityPattern);
     }
 
     double MultiscalePartitioningEngineView::getJacobianMatrixEntry(
@@ -811,6 +842,26 @@ namespace gridfire {
         return equilibrateNetwork(primingReport.primedComposition, T9, rho);
     }
 
+    bool MultiscalePartitioningEngineView::involvesSpecies(
+        const Species &species
+    ) const {
+        if (involvesSpeciesInQSE(species)) return true; // check this first since the vector is likely to be smaller so short circuit cost is less on average
+        if (involvesSpeciesInDynamic(species)) return true;
+        return false;
+    }
+
+    bool MultiscalePartitioningEngineView::involvesSpeciesInQSE(
+        const Species &species
+    ) const {
+        return std::ranges::find(m_algebraic_species, species) != m_algebraic_species.end();
+    }
+
+    bool MultiscalePartitioningEngineView::involvesSpeciesInDynamic(
+        const Species &species
+    ) const {
+        return std::ranges::find(m_dynamic_species, species) != m_dynamic_species.end();
+    }
+
     size_t MultiscalePartitioningEngineView::getSpeciesIndex(const Species &species) const {
         return m_baseEngine.getSpeciesIndex(species);
     }
@@ -1083,7 +1134,7 @@ namespace gridfire {
 
             }
 
-            if (bool group_is_coupled = (coupling_flux / leakage_flux) > FLUX_RATIO_THRESHOLD) {
+            if (coupling_flux / leakage_flux > FLUX_RATIO_THRESHOLD) {
                 LOG_TRACE_L1(
                     m_logger,
                     "Group containing {} is in equilibrium due to high coupling flux and balanced creation and destruction rate: <coupling: leakage flux = {}, coupling flux = {}, ratio = {} (Threshold: {})>, <creation: creation flux = {}, destruction flux = {}, ratio = {} order of mag (Threshold: {} order of mag)>",