feat(reaclib): working AD system and nearly working network

a few issues remain with letting the composition evolve as new species come online

src/network/private/netgraph.cpp

@@ -1,63 +1,54 @@
 #include "netgraph.h"
 #include "atomicSpecies.h"
-#include "reaclib.h"
-#include "network.h"
-#include "species.h"
 #include "const.h"
+#include "network.h"
+#include "reaclib.h"
+#include "species.h"
 
 #include "quill/LogMacros.h"
 
-#include <set>
-#include <unordered_map>
-#include <string>
-#include <vector>
-#include <stdexcept>
-#include <string_view>
 #include <cstdint>
 #include <iostream>
+#include <set>
+#include <stdexcept>
+#include <string>
+#include <string_view>
+#include <unordered_map>
+#include <vector>
 
-#include "const.h"
+#include <boost/numeric/ublas/vector.hpp>
+#include <boost/numeric/odeint.hpp>
 
 
 namespace serif::network {
-    template <typename GeneralScalarType>
+    GraphNetwork::GraphNetwork(
-    GraphNetwork<GeneralScalarType>::GraphNetwork(
         const serif::composition::Composition &composition
     ):
      Network(REACLIB),
-     m_reactions(build_reaclib_nuclear_network(composition)),
+     m_reactions(build_reaclib_nuclear_network(composition)) {
-     m_initialComposition(composition),
-     m_currentComposition(composition),
-     m_cullingThreshold(0),
-     m_T9(0) {
         syncInternalMaps();
     }
 
-    template <typename GeneralScalarType>
+    GraphNetwork::GraphNetwork(
-    GraphNetwork<GeneralScalarType>::GraphNetwork(
         const serif::composition::Composition &composition,
         const double cullingThreshold,
         const double T9
     ):
      Network(REACLIB),
-     m_reactions(build_reaclib_nuclear_network(composition, cullingThreshold, T9)),
+     m_reactions(build_reaclib_nuclear_network(composition, cullingThreshold, T9)) {
-     m_initialComposition(composition),
-     m_currentComposition(composition),
-     m_cullingThreshold(cullingThreshold),
-     m_T9(T9) {
         syncInternalMaps();
     }
 
-    template <typename GeneralScalarType>
+    void GraphNetwork::syncInternalMaps() {
-    void GraphNetwork<GeneralScalarType>::syncInternalMaps() {
         collectNetworkSpecies();
         populateReactionIDMap();
         populateSpeciesToIndexMap();
+        reserveJacobianMatrix();
+        recordADTape();
     }
 
     // --- Network Graph Construction Methods ---
-    template <typename GeneralScalarType>
+    void GraphNetwork::collectNetworkSpecies() {
-    void GraphNetwork<GeneralScalarType>::collectNetworkSpecies() {
         m_networkSpecies.clear();
         m_networkSpeciesMap.clear();
 
@@ -85,8 +76,7 @@ namespace serif::network {
 
     }
 
-    template <typename GeneralScalarType>
+    void GraphNetwork::populateReactionIDMap() {
-    void GraphNetwork<GeneralScalarType>::populateReactionIDMap() {
         LOG_INFO(m_logger, "Populating reaction ID map for REACLIB graph network (serif::network::GraphNetwork)...");
         m_reactionIDMap.clear();
         for (const auto& reaction: m_reactions) {
@@ -102,14 +92,17 @@ namespace serif::network {
         }
     }
 
-    void GraphNetwork::buildNetworkGraph() {
+    void GraphNetwork::reserveJacobianMatrix() {
-        LOG_INFO(m_logger, "Building network graph...");
+        // The implementation of this function (and others) constrains this nuclear network to a constant temperature and density during
-        m_reactions = build_reaclib_nuclear_network(m_initialComposition, m_cullingThreshold, m_T9);
+        // each evaluation.
-        syncInternalMaps();
+        size_t numSpecies = m_networkSpecies.size();
-        LOG_INFO(m_logger, "Network graph built with {} reactions and {} species.", m_reactions.size(), m_networkSpecies.size());
+        m_jacobianMatrix.clear();
+        m_jacobianMatrix.resize(numSpecies, numSpecies, false); // Sparse matrix, no initial values
+        LOG_INFO(m_logger, "Jacobian matrix resized to {} rows and {} columns.",
+                 m_jacobianMatrix.size1(), m_jacobianMatrix.size2());
     }
-        // --- Basic Accessors and Queries ---
 
+    // --- Basic Accessors and Queries ---
     const std::vector<serif::atomic::Species>& GraphNetwork::getNetworkSpecies() const {
         // Returns a constant reference to the vector of unique species in the network.
         LOG_DEBUG(m_logger, "Providing access to network species vector. Size: {}.", m_networkSpecies.size());
@@ -137,7 +130,7 @@ namespace serif::network {
         for (const auto& reactant : reaction.reactants()) {
             auto it = m_networkSpeciesMap.find(reactant.name());
             if (it != m_networkSpeciesMap.end()) {
-                stoichiometry[it->second]--; // Copy Species by value (PERF: Future performance improvements by using pointers or references)
+                stoichiometry[it->second]--; // Copy Species by value (PERF: Future performance improvements by using pointers or references (std::reference_wrapper<const ...>) or something like that)
             } else {
                 LOG_WARNING(m_logger, "Reactant species '{}' in reaction '{}' not found in network species map during stoichiometry calculation.",
                          reactant.name(), reaction.id());
@@ -148,7 +141,7 @@ namespace serif::network {
         for (const auto& product : reaction.products()) {
             auto it = m_networkSpeciesMap.find(product.name());
             if (it != m_networkSpeciesMap.end()) {
-                stoichiometry[it->second]++; // Copy Species by value (PERF: Future performance improvements by using pointers or references)
+                stoichiometry[it->second]++; // Copy Species by value (PERF: Future performance improvements by using pointers or references (std::reference_wrapper<const ...>) or something like that)
             } else {
                 LOG_WARNING(m_logger, "Product species '{}' in reaction '{}' not found in network species map during stoichiometry calculation.",
                          product.name(), reaction.id());
@@ -214,6 +207,21 @@ namespace serif::network {
         return true; // All reactions passed the conservation check
     }
 
+    void GraphNetwork::validateComposition(const serif::composition::Composition &composition, double culling, double T9) {
+
+        // Check if the requested network has already been cached.
+        // PERF: Rebuilding this should be pretty fast but it may be a good point of optimization in the future.
+        const reaclib::REACLIBReactionSet validationReactionSet = build_reaclib_nuclear_network(composition, culling, T9);
+
+
+        // This allows for dynamic network modification while retaining caching for networks which are very similar.
+        if (validationReactionSet != m_reactions) {
+            LOG_INFO(m_logger, "Reaction set not cached. Rebuilding the reaction set for T9={} and culling={}.", T9, culling);
+            m_reactions = validationReactionSet;
+            syncInternalMaps(); // Re-sync internal maps after updating reactions. Note this will also retrace the AD tape.
+        }
+    }
+
     // --- Generate Stoichiometry Matrix ---
     void GraphNetwork::generateStoichiometryMatrix() {
         LOG_INFO(m_logger, "Generating stoichiometry matrix...");
@@ -258,100 +266,191 @@ namespace serif::network {
                  m_stoichiometryMatrix.nnz()); // Assuming nnz() exists for compressed_matrix
     }
 
-    void GraphNetwork::generateJacobianMatrix() {
+    void GraphNetwork::generateJacobianMatrix(const std::vector<double> &Y, const double T9,
-        LOG_INFO(m_logger, "Generating jacobian matrix...");
+        const double rho) {
 
-        // Task 1: Set dimensions and initialize the matrix
+        LOG_INFO(m_logger, "Generating jacobian matrix for T9={}, rho={}..", T9, rho);
-        size_t numSpecies = m_networkSpecies.size();
-        m_stoichiometryMatrix.resize(numSpecies, numSpecies, false);
-
-        LOG_INFO(m_logger, "Jacobian matrix initialized with dimensions: {} rows (species) x {} columns (species).",
-                 numSpecies, numSpecies);
-    }
-
-    std::vector<double> GraphNetwork::calculateRHS(const std::vector<double>& Y, const double T9, const double rho) const {
-        std::vector<double> dotY(m_networkSpecies.size(), 0.0);
-        const size_t numReactions = m_reactions.size();
         const size_t numSpecies = m_networkSpecies.size();
 
-        if (m_stoichiometryMatrix.size1() != numSpecies || m_stoichiometryMatrix.size2() != numReactions) {
+        // 1. Pack the input variables into a vector for CppAD
-            LOG_ERROR(m_logger, "Stoichiometry matrix dimensions do not match network species and reactions sizes.");
+        std::vector<double> adInput(numSpecies + 2, 0.0); // +2 for T9 and rho
-            throw std::runtime_error("Stoichiometry matrix dimensions mismatch.");
+        for (size_t i = 0; i < numSpecies; ++i) {
+            adInput[i] = Y[i];
         }
+        adInput[numSpecies]     = T9;  // T9
+        adInput[numSpecies + 1] = rho; // rho
 
-        for (size_t reactionIndex = 0; reactionIndex < numReactions; ++reactionIndex) {
+        // 2. Calculate the full jacobian
-            const auto& currentReaction = m_reactions[reactionIndex];
+        const std::vector<double> dotY = m_rhsADFun.Jacobian(adInput);
-            const double reactionRatePerVolumePerTime = calculateReactionRate(currentReaction, Y, T9, rho);
 
-            // dY_i/dt = ∑ υ_ij * R_j * mass_in_grams
+        // 3. Pack jacobian vector into sparse matrix
-            // TODO: make sure the unit conversion is correct after calculate reaction rate.
+        m_jacobianMatrix.clear();
-
+        for (size_t i = 0; i < numSpecies; ++i) {
-            for (size_t speciesIndex = 0; speciesIndex < numSpecies; ++speciesIndex) {
+            for (size_t j = 0; j < numSpecies; ++j) {
-                const int nu_ij = m_stoichiometryMatrix(speciesIndex, reactionIndex);
+                const double value = dotY[i * (numSpecies + 2) + j];
-                if (nu_ij != 0) {
+                if (std::abs(value) > MIN_JACOBIAN_THRESHOLD) {
-                    const double speciesAtomicMassAMU = m_networkSpecies[speciesIndex].mass();
+                    m_jacobianMatrix(i, j) = value;
-                    dotY[speciesIndex] += (nu_ij * reactionRatePerVolumePerTime * speciesAtomicMassAMU)/rho;
                 }
             }
         }
-        return dotY;
+        LOG_INFO(m_logger, "Jacobian matrix generated with dimensions: {} rows x {} columns.", m_jacobianMatrix.size1(), m_jacobianMatrix.size2());
     }
 
-    double GraphNetwork::calculateReactionRate(const reaclib::REACLIBReaction &reaction, const std::vector<double> &Y,
+    void GraphNetwork::detectStiff(const NetIn &netIn, const double T9, const double numSpecies, const boost::numeric::ublas::vector<double>& Y) {
-        const double T9, const double rho) const {
+        // --- Heuristic for automatic stiffness detection ---
-        const auto &constants  = serif::constant::Constants::getInstance();
+        const std::vector<double> initial_y_stl(Y.begin(), Y.begin() + numSpecies); // Copy only the species abundances, not the specific energy rate
+        const auto derivatives = calculateAllDerivatives<double>(initial_y_stl, T9, netIn.density);
+        const std::vector<double>& initial_dotY = derivatives.dydt;
 
-        const auto u = constants.get("u"); // Atomic mass unit in g/mol
+        double min_timescale = std::numeric_limits<double>::max();
-        const double k_reaction = reaction.calculate_rate(T9); // PERF: Consider precomputing all of these and putting them into an O(1) lookup table.
+        double max_timescale = 0.0;
 
-        double reactant_product = 1.0;
+        for (size_t i = 0; i < numSpecies; ++i) {
-
+            if (std::abs(initial_dotY[i]) > 0) {
-        std::unordered_map<std::string, int> reactant_counts;
+                const double timescale = std::abs(Y(i) / initial_dotY[i]);
-        reactant_counts.reserve(reaction.reactants().size());
+                if (timescale > max_timescale) {max_timescale = timescale;}
-        for (const auto& reactant : reaction.reactants()) {
+                if (timescale < min_timescale) {min_timescale = timescale;}
-            reactant_counts[std::string(reactant.name())]++;
+            }
         }
 
-        for (const auto& [species_name, count] : reactant_counts) {
+        const double stiffnessRatio = max_timescale / min_timescale;
-            constexpr double minThreshold = 1e-18;
-            auto species_it = m_speciesToIndexMap.find(m_networkSpeciesMap.at(species_name));
-            if (species_it == m_speciesToIndexMap.end()) {
-                LOG_ERROR(m_logger, "Reactant species '{}' not found in species to index map for reaction '{}'.",
-                         species_name, reaction.id());
-                throw std::runtime_error("Reactant species not found in species to index map: " + species_name);
-            }
-            const size_t species_index = species_it->second;
-            const double Yi = Y[species_index];
-            double Ai = m_networkSpecies[species_index].a();
 
-            if (Yi < minThreshold) {
+        // TODO: Pull this out into a configuration option or a more sophisticated heuristic.
-                return 0.0; // If any reactant is below a threshold, return zero rate.
+        constexpr double stiffnessThreshold = 1.0e6; // This is a heuristic threshold, can be tuned based on network characteristics.
-            }
 
-            double atomicMassAMU = m_networkSpecies[species_index].mass();
+        LOG_INFO(m_logger, "Stiffness ratio is {} (max timescale: {}, min timescale: {}).", stiffnessRatio, max_timescale, min_timescale);
-
+        if (stiffnessRatio > stiffnessThreshold) {
-            // Convert to number density
+            LOG_INFO(m_logger, "Network is detected to be stiff. Using stiff ODE solver.");
-            double ni;
+            m_stiff = true;
-            const double denominator = atomicMassAMU * u.value;
-            if (denominator > minThreshold) {
-                ni = (Yi * rho) / (denominator);
         } else {
-                ni = 0.0;
+            LOG_INFO(m_logger, "Network is detected to be non-stiff. Using non-stiff ODE solver.");
+            m_stiff = false;
         }
-
-            reactant_product *= ni;
-
-            // Apply factorial correction for identical reactions
-            if (count > 1) {
-                reactant_product /= static_cast<double>(std::tgamma(count + 1)); // Gamma function for factorial
-            }
-        }
-        const double Na = constants.get("N_a").value; // Avogadro's number in mol^-1
-        const double molarCorrectionFactor = std::pow(Na, reaction.reactants().size() - 1);
-        return (reactant_product * k_reaction) / molarCorrectionFactor; // reaction rate in per volume per time (particles/cm^3/s)
     }
 
     NetOut GraphNetwork::evaluate(const NetIn &netIn) {
-        return Network::evaluate(netIn);
+        namespace ublas = boost::numeric::ublas;
+        namespace odeint = boost::numeric::odeint;
+
+        const double T9 = netIn.temperature / 1e9; // Convert temperature from Kelvin to T9 (T9 = T / 1e9)
+        validateComposition(netIn.composition, netIn.culling, T9);
+
+        const double numSpecies = m_networkSpecies.size();
+        constexpr double abs_tol = 1.0e-8;
+        constexpr double rel_tol = 1.0e-8;
+
+        int stepCount = 0;
+
+        // TODO: Pull these out into configuration options
+
+        ODETerm rhs_functor(*this, T9, netIn.density);
+
+
+        ublas::vector<double> Y(numSpecies + 1);
+        for (size_t i = 0; i < numSpecies; ++i) {
+            const auto& species = m_networkSpecies[i];
+            // Get the mass fraction for this specific species from the input object
+            Y(i) = netIn.composition.getMassFraction(std::string(species.name()));
+        }
+        Y(numSpecies) = 0; // initial specific energy rate, will be updated later
+
+        detectStiff(netIn, T9, numSpecies, Y);
+
+        if (m_stiff) {
+            JacobianTerm jacobian_functor(*this, T9, netIn.density);
+            LOG_INFO(m_logger, "Making use of stiff ODE solver for network evaluation.");
+            auto stepper = odeint::make_controlled<odeint::rosenbrock4<double>>(abs_tol, rel_tol);
+            stepCount = odeint::integrate_adaptive(
+                stepper,
+                std::make_pair(rhs_functor, jacobian_functor),
+                Y,
+                0.0, // Start time
+                netIn.tMax,
+                netIn.dt0
+                );
+
+        } else {
+            LOG_INFO(m_logger, "Making use of ODE solver (non-stiff) for network evaluation.");
+            using state_type = ublas::vector<double>;
+            auto stepper = odeint::make_controlled<odeint::runge_kutta_dopri5<state_type>>(abs_tol, rel_tol);
+            stepCount = odeint::integrate_adaptive(
+                stepper,
+                rhs_functor,
+                Y,
+                0.0, // Start time
+                netIn.tMax,
+                netIn.dt0
+            );
+
+
+        }
+
+        double sumY = 0.0;
+        for (int i = 0; i < numSpecies; ++i) { sumY += Y(i); }
+        for (int i = 0; i < numSpecies; ++i) { Y(i) /= sumY; }
+
+        // --- Marshall output variables ---
+        // PERF: Im sure this step could be tuned to avoid so many copies, that is a job for another day
+        std::vector<std::string> speciesNames;
+        speciesNames.reserve(numSpecies);
+        for (const auto& species : m_networkSpecies) {
+            speciesNames.push_back(std::string(species.name()));
+        }
+
+        std::vector<double> finalAbundances(Y.begin(), Y.begin() + numSpecies);
+        serif::composition::Composition outputComposition(speciesNames, finalAbundances);
+        outputComposition.finalize(true);
+
+        NetOut netOut;
+        netOut.composition = outputComposition;
+        netOut.num_steps = stepCount;
+        netOut.energy = Y(numSpecies); // The last element in Y is the specific energy rate
+
+        return netOut;
+
+    }
+
+    void GraphNetwork::recordADTape() {
+        LOG_INFO(m_logger, "Recording AD tape for the RHS calculation...");
+
+        // Task 1: Set dimensions and initialize the matrix
+        const size_t numSpecies = m_networkSpecies.size();
+        if (numSpecies == 0) {
+            LOG_ERROR(m_logger, "Cannot record AD tape: No species in the network.");
+            throw std::runtime_error("Cannot record AD tape: No species in the network.");
+        }
+        const size_t numADInputs = numSpecies + 2; // Note here that by not letting T9 and rho be independent variables, we are constraining the network to a constant temperature and density during each evaluation.
+
+        // --- CppAD Tape Recording ---
+        // 1. Declare independent variable (adY)
+        //    We also initialize the dummy variable for tape recording (these tell CppAD what the derivative chain looks like).
+        //    Their numeric values are irrelevant except for in so far as they avoid numerical instabilities.
+
+        // Distribute total mass fraction uniformly between species in the dummy variable space
+        const auto uniformMassFraction = static_cast<CppAD::AD<double>>(1.0 / numSpecies);
+        std::vector<CppAD::AD<double>> adInput(numADInputs, uniformMassFraction);
+        adInput[numSpecies]     = 1.0; // Dummy T9
+        adInput[numSpecies + 1] = 1.0; // Dummy rho
+
+        // 3. Declare independent variables (what CppAD will differentiate wrt.)
+        //    This also beings the tape recording process.
+        CppAD::Independent(adInput);
+
+        std::vector<CppAD::AD<double>> adY(numSpecies);
+        for(size_t i = 0; i < numSpecies; ++i) {
+            adY[i] = adInput[i];
+        }
+        const CppAD::AD<double> adT9  = adInput[numSpecies];
+        const CppAD::AD<double> adRho = adInput[numSpecies + 1];
+
+
+        // 5. Call the actual templated function
+        // We let T9 and rho be constant, so we pass them as fixed values.
+        auto derivatives = calculateAllDerivatives<CppAD::AD<double>>(adY, adT9, adRho);
+
+        m_rhsADFun.Dependent(adInput, derivatives.dydt);
+
+        LOG_INFO(m_logger, "AD tape recorded successfully for the RHS calculation. Number of independent variables: {}.",
+                 adInput.size());
     }
 }

src/network/private/network.cpp

@@ -33,7 +33,8 @@ namespace serif::network {
         m_config(serif::config::Config::getInstance()),
         m_logManager(serif::probe::LogManager::getInstance()),
         m_logger(m_logManager.getLogger("log")),
-        m_format(format) {
+        m_format(format),
+        m_constants(serif::constant::Constants::getInstance()){
         if (format == NetworkFormat::UNKNOWN) {
             LOG_ERROR(m_logger, "nuclearNetwork::Network::Network() called with UNKNOWN format");
             throw std::runtime_error("nuclearNetwork::Network::Network() called with UNKNOWN format");
@@ -73,6 +74,7 @@ namespace serif::network {
     serif::network::reaclib::REACLIBReactionSet build_reaclib_nuclear_network(const serif::composition::Composition &composition) {
         using namespace serif::network::reaclib;
         REACLIBReactionSet reactions;
+        auto logger = serif::probe::LogManager::getInstance().getLogger("log");
 
         if (!s_initialized) {
             LOG_INFO(serif::probe::LogManager::getInstance().getLogger("log"), "REACLIB reactions not initialized. Calling initializeAllReaclibReactions()...");
@@ -89,6 +91,7 @@ namespace serif::network {
                 }
             }
             if (gotReaction) {
+                LOG_INFO(logger, "Adding reaction {} to REACLIB reaction set.", reaction.id());
                 reactions.add_reaction(reaction);
             }
         }

src/network/public/approx8.h

@@ -306,13 +306,13 @@ namespace serif::network::approx8{
          * @brief Sets whether the solver should use a stiff method.
          * @param stiff Boolean indicating if a stiff method should be used.
          */
-        void setStiff(bool stiff);
+        void setStiff(bool stiff) override;
 
         /**
          * @brief Checks if the solver is using a stiff method.
          * @return Boolean indicating if a stiff method is being used.
          */
-        bool isStiff() const { return m_stiff; }
+        bool isStiff() const override { return m_stiff; }
     private:
         vector_type m_y;
         double m_tMax = 0;

src/network/public/netgraph.h

@@ -1,16 +1,20 @@
 #pragma once
 
-#include "network.h"
-#include "reaclib.h"
 #include "atomicSpecies.h"
 #include "composition.h"
+#include "network.h"
+#include "reaclib.h"
 
+#include <string>
 #include <unordered_map>
 #include <vector>
-#include <string>
 
 #include <boost/numeric/ublas/matrix_sparse.hpp>
 
+#include "cppad/cppad.hpp"
+
+#include "quill/LogMacros.h"
+
 
 // PERF: The function getNetReactionStoichiometry returns a map of species to their stoichiometric coefficients for a given reaction.
 //       this makes extra copies of the species, which is not ideal and could be optimized further.
@@ -18,54 +22,248 @@
 //       REACLIBReactions are quite large data structures, so this could be a performance bottleneck.
 
 namespace serif::network {
-    template <typename GeneralScalarType>
+    // Define a concept to check if a type is one of our allowed scalar types
+    template<typename T>
+    concept IsArithmeticOrAD = std::is_same_v<T, double> || std::is_same_v<T, CppAD::AD<double>>;
+
+    static constexpr double MIN_DENSITY_THRESHOLD = 1e-18;
+    static constexpr double MIN_ABUNDANCE_THRESHOLD = 1e-18;
+    static constexpr double MIN_JACOBIAN_THRESHOLD = 1e-24;
+
     class GraphNetwork final : public Network {
     public:
         explicit GraphNetwork(const serif::composition::Composition &composition);
-        explicit GraphNetwork(const serif::composition::Composition &composition, const double cullingThreshold, const double T9);
+        explicit GraphNetwork(const serif::composition::Composition &composition,
+                              const double cullingThreshold, const double T9);
 
         NetOut evaluate(const NetIn &netIn) override;
 
         [[nodiscard]] const std::vector<serif::atomic::Species>& getNetworkSpecies() const;
         [[nodiscard]] const reaclib::REACLIBReactionSet& getNetworkReactions() const;
-        [[nodiscard]] std::unordered_map<serif::atomic::Species, int> getNetReactionStoichiometry(const reaclib::REACLIBReaction& reaction) const;
+        [[nodiscard]] std::unordered_map<serif::atomic::Species, int> getNetReactionStoichiometry(
+            const reaclib::REACLIBReaction &reaction) const;
 
         [[nodiscard]] bool involvesSpecies(const serif::atomic::Species& species) const;
 
-        std::vector<std::vector<std::string>> detectCycles() const;
+        [[deprecated("not implimented")]] std::vector<std::vector<std::string>> detectCycles() const;
-
+        [[deprecated("not implimented")]] std::vector<std::string> topologicalSortReactions() const;
-        std::vector<std::string> topologicalSortReactions() const;
     private:
         reaclib::REACLIBReactionSet m_reactions; ///< Set of REACLIB reactions for the network.
-        serif::composition::Composition m_initialComposition;
-        serif::composition::Composition m_currentComposition;
-        GeneralScalarType m_cullingThreshold; ///< Threshold for culling reactions.
-        GeneralScalarType m_T9; ///< Temperature in T9 units to evaluate culling.
 
         std::vector<serif::atomic::Species> m_networkSpecies; ///< The species in the network.
         std::unordered_map<std::string_view, serif::atomic::Species> m_networkSpeciesMap;
         std::unordered_map<std::string_view, const reaclib::REACLIBReaction> m_reactionIDMap; ///< Map of reaction IDs to REACLIB reactions.
 
         boost::numeric::ublas::compressed_matrix<int> m_stoichiometryMatrix; ///< Stoichiometry matrix for the network.
-        boost::numeric::ublas::compressed_matrix<GeneralScalarType> m_jacobianMatrix; ///< Jacobian matrix for the network.
+        boost::numeric::ublas::compressed_matrix<double> m_jacobianMatrix; ///< Jacobian matrix for the network.
         std::unordered_map<serif::atomic::Species, size_t> m_speciesToIndexMap; ///< Map of species to their index in the stoichiometry matrix.
 
+        CppAD::ADFun<double> m_rhsADFun; ///< AD tape
+
+        struct ODETerm {
+            GraphNetwork& m_network; ///< Reference to the network
+            const double m_T9;
+            const double m_rho;
+            const size_t m_numSpecies;
+
+            ODETerm(GraphNetwork& network, const double T9, double rho) :
+             m_network(network), m_T9(T9), m_rho(rho), m_numSpecies(m_network.m_networkSpecies.size()) {}
+
+            void operator()(const boost::numeric::ublas::vector<double>&Y, boost::numeric::ublas::vector<double>& dYdt, double /* t */) {
+                const std::vector<double> y(Y.begin(), Y.begin() + m_numSpecies);
+
+                StepDerivatives<double> derivatives = m_network.calculateAllDerivatives<double>(y, m_T9, m_rho);
+
+                dYdt.resize(m_numSpecies + 1);
+                std::copy(derivatives.dydt.begin(), derivatives.dydt.end(), dYdt.begin());
+                dYdt(m_numSpecies) = derivatives.specificEnergyRate; // Last element is the specific energy rate
+            }
+        };
+
+        struct JacobianTerm {
+            GraphNetwork& m_network;
+            const double m_T9;
+            const double m_rho;
+            const size_t m_numSpecies;
+
+            JacobianTerm(GraphNetwork& network, const double T9, const double rho) :
+                m_network(network), m_T9(T9), m_rho(rho), m_numSpecies(m_network.m_networkSpecies.size()) {}
+
+            void operator()(const boost::numeric::ublas::vector<double>& Y, boost::numeric::ublas::matrix<double>& J, double /* t */, boost::numeric::ublas::vector<double>& /* dfdt */) {
+                const std::vector<double> y_species(Y.begin(), Y.begin() + m_numSpecies);
+
+                m_network.generateJacobianMatrix(y_species, m_T9, m_rho);
+
+                J.resize(m_numSpecies + 1, m_numSpecies + 1);
+                J.clear(); // Zero out the matrix
+
+                // PERF: Could probably be optimized
+                for (auto it1 = m_network.m_jacobianMatrix.begin1(); it1 != m_network.m_jacobianMatrix.end1(); ++it1) {
+                    for (auto it2 = it1.begin(); it2 != it1.end(); ++it2) {
+                        J(it2.index1(), it2.index2()) = *it2;
+                    }
+                }
+            }
+        };
+
+        template <IsArithmeticOrAD T>
+        struct StepDerivatives {
+            std::vector<T> dydt;
+            T specificEnergyRate = T(0.0);
+        };
+
 
     private:
         void syncInternalMaps();
         void collectNetworkSpecies();
         void populateReactionIDMap();
         void populateSpeciesToIndexMap();
-        void buildNetworkGraph();
+        void reserveJacobianMatrix();
+        void recordADTape();
+
+        // --- Validation methods ---
         bool validateConservation() const;
 
+        void validateComposition(const serif::composition::Composition &composition, double culling, double T9);
+
+        // --- Simple Derivative Calculations ---
+        template <IsArithmeticOrAD T>
+        StepDerivatives<T> calculateAllDerivatives(const std::vector<T>& Y, T T9, T rho) const;
+
         // --- Generate the system matrices ---
         void generateStoichiometryMatrix();
-        void generateJacobianMatrix();
+        void generateJacobianMatrix(const std::vector<double>& Y, const double T9, const double rho);
 
-        std::vector<GeneralScalarType> calculateRHS(const std::vector<GeneralScalarType>& Y, const GeneralScalarType T9, const GeneralScalarType rho) const;
+        template <IsArithmeticOrAD GeneralScalarType>
-        GeneralScalarType calculateReactionRate(const reaclib::REACLIBReaction &reaction, const std::vector<GeneralScalarType>& Y, const GeneralScalarType T9, const GeneralScalarType rho) const;
+        std::vector<GeneralScalarType> calculateRHS(const std::vector<GeneralScalarType> &Y, const GeneralScalarType T9,
+                                                    const GeneralScalarType rho) const;
+        template <IsArithmeticOrAD GeneralScalarType>
+        GeneralScalarType calculateReactionRate(const reaclib::REACLIBReaction &reaction,
+                                                const std::vector<GeneralScalarType> &Y, const GeneralScalarType T9,
+                                                const GeneralScalarType rho) const;
 
-        void pruneNetworkGraph();
+        void detectStiff(const NetIn &netIn, double T9, double numSpecies, const boost::numeric::ublas::vector<double>& Y);
     };
+
+
+    template<IsArithmeticOrAD T>
+    typename GraphNetwork::template StepDerivatives<T> GraphNetwork::calculateAllDerivatives(
+        const std::vector<T> &Y, T T9, T rho) const {
+        StepDerivatives<T> result;
+        result.dydt.resize(m_networkSpecies.size(), static_cast<T>(0.0));
+
+        if (rho < MIN_DENSITY_THRESHOLD) {
+            return result; // Return zero for everything if density is too low
+        }
+
+        const T u = static_cast<T>(m_constants.get("u").value); // Atomic mass unit in grams
+        const T MeV_to_erg = static_cast<T>(m_constants.get("MeV_to_erg").value);
+
+        T volumetricEnergyRate = static_cast<T>(0.0); // Accumulator for erg / (cm^3 * s)
+
+        // --- SINGLE LOOP OVER ALL REACTIONS ---
+        for (size_t reactionIndex = 0; reactionIndex < m_reactions.size(); ++reactionIndex) {
+            const auto& reaction = m_reactions[reactionIndex];
+
+            // 1. Calculate reaction rate
+            const T reactionRate = calculateReactionRate(reaction, Y, T9, rho);
+
+            // 2. Use the rate to update all relevant species derivatives (dY/dt)
+            for (size_t speciesIndex = 0; speciesIndex < m_networkSpecies.size(); ++speciesIndex) {
+                const T nu_ij = static_cast<T>(m_stoichiometryMatrix(speciesIndex, reactionIndex));
+                if (nu_ij != 0) {
+                    const T speciesAtomicMassAMU = static_cast<T>(m_networkSpecies[speciesIndex].mass());
+                    const T speciesAtomicMassGrams = speciesAtomicMassAMU * u;
+                    result.dydt[speciesIndex] += (nu_ij * reactionRate * speciesAtomicMassGrams) / rho;
+                }
+            }
+
+            // 3. Use the same rate to update the energy generation rate
+            const T q_value_ergs = static_cast<T>(reaction.qValue()) * MeV_to_erg;
+            volumetricEnergyRate += reactionRate * q_value_ergs;
+        }
+
+        result.specificEnergyRate = volumetricEnergyRate / rho;
+
+        return result;
+    }
+
+
+    template <IsArithmeticOrAD GeneralScalarType>
+    GeneralScalarType GraphNetwork::calculateReactionRate(const reaclib::REACLIBReaction &reaction, const std::vector<GeneralScalarType> &Y,
+        const GeneralScalarType T9, const GeneralScalarType rho) const {
+        const auto &constants  = serif::constant::Constants::getInstance();
+
+        const auto u = constants.get("u"); // Atomic mass unit in g/mol
+        const GeneralScalarType uValue = static_cast<GeneralScalarType>(u.value); // Convert to double for calculations
+        const GeneralScalarType k_reaction = reaction.calculate_rate(T9);
+
+        GeneralScalarType reactant_product_or_particle_densities = static_cast<GeneralScalarType>(1.0);
+
+        std::unordered_map<std::string, int> reactant_counts;
+        reactant_counts.reserve(reaction.reactants().size());
+        for (const auto& reactant : reaction.reactants()) {
+            reactant_counts[std::string(reactant.name())]++;
+        }
+
+        const GeneralScalarType minAbundanceThreshold = static_cast<GeneralScalarType>(MIN_ABUNDANCE_THRESHOLD);
+        const GeneralScalarType minDensityThreshold = static_cast<GeneralScalarType>(MIN_DENSITY_THRESHOLD);
+
+        if (rho < minDensityThreshold) {
+            // LOG_INFO(m_logger, "Density is below the minimum threshold ({} g/cm^3), returning zero reaction rate for reaction '{}'.",
+            //          CppAD::Value(rho), reaction.id()); // CppAD::Value causes a compilation error here, not clear why...
+            return static_cast<GeneralScalarType>(0.0); // If density is below a threshold, return zero rate.
+        }
+
+        for (const auto& [species_name, count] : reactant_counts) {
+            auto species_it = m_speciesToIndexMap.find(m_networkSpeciesMap.at(species_name));
+            if (species_it == m_speciesToIndexMap.end()) {
+                LOG_ERROR(m_logger, "Reactant species '{}' not found in species to index map for reaction '{}'.",
+                         species_name, reaction.id());
+                throw std::runtime_error("Reactant species not found in species to index map: " + species_name);
+            }
+            const size_t species_index = species_it->second;
+            const GeneralScalarType Yi = Y[species_index];
+
+            if (Yi < minAbundanceThreshold) {
+                return static_cast<GeneralScalarType>(0.0); // If any reactant is below a threshold, return zero rate.
+            }
+
+            GeneralScalarType atomicMassAMU = static_cast<GeneralScalarType>(m_networkSpecies[species_index].mass());
+
+            // Convert to number density
+            GeneralScalarType ni;
+            const GeneralScalarType denominator = atomicMassAMU * uValue;
+            if (denominator > minDensityThreshold) {
+                ni = (Yi * rho) / (denominator);
+            } else {
+                ni = static_cast<GeneralScalarType>(0.0);
+            }
+
+            reactant_product_or_particle_densities *= ni;
+
+            // Apply factorial correction for identical reactions
+            if (count > 1) {
+                reactant_product_or_particle_densities /= static_cast<GeneralScalarType>(std::tgamma(static_cast<double>(count + 1))); // Gamma function for factorial
+            }
+        }
+        const GeneralScalarType Na = static_cast<GeneralScalarType>(constants.get("N_a").value); // Avogadro's number in mol^-1
+        const int numReactants = static_cast<int>(reaction.reactants().size());
+        GeneralScalarType molarCorrectionFactor = static_cast<GeneralScalarType>(1.0); // No correction needed for single reactant reactions
+        if (numReactants > 1) {
+            molarCorrectionFactor = CppAD::pow(Na, static_cast<GeneralScalarType>(reaction.reactants().size() - 1));
+        }
+        return (reactant_product_or_particle_densities * k_reaction) / molarCorrectionFactor; // reaction rate in per volume per time (particles/cm^3/s)
+    }
+
+    template <IsArithmeticOrAD GeneralScalarType>
+    std::vector<GeneralScalarType> GraphNetwork::calculateRHS(
+        const std::vector<GeneralScalarType>& Y,
+        const GeneralScalarType T9,
+        const GeneralScalarType rho) const {
+
+        auto derivatives = calculateAllDerivatives<GeneralScalarType>(Y, T9, rho);
+        return derivatives.dydt;
+    }
+
 };

src/network/public/network.h

@@ -29,6 +29,8 @@
 #include "reaclib.h"
 #include <unordered_map>
 
+#include "const.h"
+
 namespace serif::network {
 
     enum NetworkFormat {
@@ -67,6 +69,7 @@ namespace serif::network {
         double temperature; ///< Temperature in Kelvin
         double density; ///< Density in g/cm^3
         double energy; ///< Energy in ergs
+        double culling = 0.0; ///< Culling threshold for reactions (default is 0.0, meaning no culling)
     };
 
     /**
@@ -87,6 +90,11 @@ namespace serif::network {
         serif::composition::Composition composition; ///< Composition of the network after evaluation
         int num_steps; ///< Number of steps taken in the evaluation
         double energy; ///< Energy in ergs after evaluation
+
+        friend std::ostream& operator<<(std::ostream& os, const NetOut& netOut) {
+            os << "NetOut(composition=" << netOut.composition << ", num_steps=" << netOut.num_steps << ", energy=" << netOut.energy << ")";
+            return os;
+        }
     };
 
     /**
@@ -119,12 +127,18 @@ namespace serif::network {
              */
             virtual NetOut evaluate(const NetIn &netIn);
 
+            virtual bool isStiff() const { return m_stiff; }
+            virtual void setStiff(const bool stiff) { m_stiff = stiff; }
+
         protected:
             serif::config::Config& m_config; ///< Configuration instance
             serif::probe::LogManager& m_logManager; ///< Log manager instance
             quill::Logger* m_logger; ///< Logger instance
 
             NetworkFormat m_format; ///< Format of the network
+            serif::constant::Constants& m_constants;
+
+            bool m_stiff = false; ///< Flag indicating if the network is stiff
     };
 
     class ReaclibNetwork final : public Network {

tests/network/approx8Test.cpp

@@ -6,7 +6,7 @@
 #include "network.h"
 #include "composition.h"
 #include "reaclib.h"
-#include "reactions.h"
+#include "netgraph.h"
 
 #include <vector>
 
@@ -74,8 +74,19 @@ TEST_F(approx8Test, reaclib) {
     serif::composition::Composition composition;
     composition.registerSymbol(symbols, true);
     composition.setMassFraction(symbols, comp);
-    bool isFinalized = composition.finalize(true);
+    composition.finalize(true);
 
-    reaclib::REACLIBReactionSet reactionSet = build_reaclib_nuclear_network(composition);
+    NetIn netIn;
-    std::cout << reactionSet.size() << std::endl;
+    netIn.composition = composition;
+    netIn.temperature = 1e7;
+    netIn.density = 1e2;
+    netIn.energy = 0.0;
+
+    netIn.tMax = 3.15e17;
+    netIn.dt0 = 1e12;
+
+    GraphNetwork network(composition);
+    NetOut netOut;
+    netOut = network.evaluate(netIn);
+    std::cout << netOut << std::endl;
 }