diff --git a/Doxyfile b/Doxyfile
index bb33f630..5cd08c0b 100644
--- a/Doxyfile
+++ b/Doxyfile
@@ -48,7 +48,7 @@ PROJECT_NAME           = GridFire
 # could be handy for archiving the generated documentation or if some version
 # control system is used.
 
-PROJECT_NUMBER         = v0.7.0-alpha
+PROJECT_NUMBER         = v0.7.0_alpha_2025_10_25
 
 # Using the PROJECT_BRIEF tag one can provide an optional one line description
 # for a project that appears at the top of each page and should give viewers a
diff --git a/build-python/meson.build b/build-python/meson.build
index 671ed52c..104ff5f1 100644
--- a/build-python/meson.build
+++ b/build-python/meson.build
@@ -1,5 +1,5 @@
 # --- Python Extension Setup ---
-py_installation = import('python').find_installation('python3')
+py_installation = import('python').find_installation('python3', pure: false)
 
 gridfire_py_deps = [
     pybind11_dep,
@@ -10,11 +10,10 @@ gridfire_py_deps = [
 ]
 
 py_mod = py_installation.extension_module(
-    'gridfire', # Name of the generated .so/.pyd file (without extension)
+    '_gridfire', # Name of the generated .so/.pyd file (without extension)
     sources: [
         meson.project_source_root() + '/src/python/bindings.cpp',
         meson.project_source_root() + '/src/python/types/bindings.cpp',
-        meson.project_source_root() + '/src/python/expectations/bindings.cpp',
         meson.project_source_root() + '/src/python/partition/bindings.cpp',
         meson.project_source_root() + '/src/python/partition/trampoline/py_partition.cpp',
         meson.project_source_root() + '/src/python/reaction/bindings.cpp',
@@ -27,9 +26,52 @@ py_mod = py_installation.extension_module(
         meson.project_source_root() + '/src/python/engine/trampoline/py_engine.cpp',
         meson.project_source_root() + '/src/python/solver/bindings.cpp',
         meson.project_source_root() + '/src/python/solver/trampoline/py_solver.cpp',
+        meson.project_source_root() + '/src/python/policy/bindings.cpp',
+        meson.project_source_root() + '/src/python/policy/trampoline/py_policy.cpp',
         meson.project_source_root() + '/src/python/utils/bindings.cpp',
     ],
     dependencies : gridfire_py_deps,
     cpp_args : ['-UNDEBUG'], # Example: Ensure assertions are enabled if needed
     install : true,
-)
\ No newline at end of file
+    subdir: 'gridfire',
+)
+
+
+py_installation.install_sources(
+    files(
+        meson.project_source_root() + '/src/python/gridfire/__init__.py',
+        meson.project_source_root() + '/stubs/gridfire/_gridfire/__init__.pyi',
+        meson.project_source_root() + '/stubs/gridfire/_gridfire/exceptions.pyi',
+        meson.project_source_root() + '/stubs/gridfire/_gridfire/partition.pyi',
+        meson.project_source_root() + '/stubs/gridfire/_gridfire/reaction.pyi',
+        meson.project_source_root() + '/stubs/gridfire/_gridfire/screening.pyi',
+        meson.project_source_root() + '/stubs/gridfire/_gridfire/io.pyi',
+        meson.project_source_root() + '/stubs/gridfire/_gridfire/solver.pyi',
+        meson.project_source_root() + '/stubs/gridfire/_gridfire/policy.pyi',
+        meson.project_source_root() + '/stubs/gridfire/_gridfire/type.pyi'
+    ),
+    subdir: 'gridfire',
+)
+
+py_installation.install_sources(
+    files(
+        meson.project_source_root() + '/stubs/gridfire/_gridfire/engine/__init__.pyi',
+        meson.project_source_root() + '/stubs/gridfire/_gridfire/engine/diagnostics.pyi',
+    ),
+    subdir: 'gridfire/engine',
+)
+
+py_installation.install_sources(
+    files(
+        meson.project_source_root() + '/stubs/gridfire/_gridfire/utils/__init__.pyi',
+    ),
+    subdir: 'gridfire/utils',
+)
+
+py_installation.install_sources(
+    files(
+        meson.project_source_root() + '/stubs/gridfire/_gridfire/utils/hashing/__init__.pyi',
+        meson.project_source_root() + '/stubs/gridfire/_gridfire/utils/hashing/reaction.pyi',
+    ),
+    subdir: 'gridfire/utils/hashing',
+)
diff --git a/meson.build b/meson.build
index 80ca4e02..e41e25fa 100644
--- a/meson.build
+++ b/meson.build
@@ -18,7 +18,7 @@
 #   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 #
 # *********************************************************************** #
-project('GridFire', 'cpp', version: 'v0.7.0-alpha', default_options: ['cpp_std=c++23'], meson_version: '>=1.5.0')
+project('GridFire', 'cpp', version: 'v0.7.0_alpha_2025_10_25', default_options: ['cpp_std=c++23'], meson_version: '>=1.5.0')
 
 # Add default visibility for all C++ targets
 add_project_arguments('-fvisibility=default', language: 'cpp')
diff --git a/pip_install_mac_patch.sh b/pip_install_mac_patch.sh
index 4f5b7508..fd0747b0 100755
--- a/pip_install_mac_patch.sh
+++ b/pip_install_mac_patch.sh
@@ -67,10 +67,10 @@ echo ""
 echo -e "${GREEN}Step 3: Locating installed gridfire extension module...${NC}"
 
 # Find the .so file
-SO_FILE=$(find "$SITE_PACKAGES" -name "gridfire.cpython-*-darwin.so" 2>/dev/null | head -n 1)
+SO_FILE=$(find "$SITE_PACKAGES/gridfire" -name "_gridfire.cpython-*-darwin.so" 2>/dev/null | head -n 1)
 
 if [ -z "$SO_FILE" ]; then
-    echo -e "${RED}Error: Could not find gridfire.cpython-*-darwin.so in $SITE_PACKAGES/gridfire${NC}"
+    echo -e "${RED}Error: Could not find _gridfire.cpython-*-darwin.so in $SITE_PACKAGES/gridfire${NC}"
     echo "Installation may have failed or the file is in an unexpected location."
     exit 1
 fi
diff --git a/src/include/gridfire/engine/engine.h b/src/include/gridfire/engine/engine.h
index 028315a9..3be40c1c 100644
--- a/src/include/gridfire/engine/engine.h
+++ b/src/include/gridfire/engine/engine.h
@@ -179,4 +179,5 @@
 
 #include "gridfire/engine/views/engine_views.h"
 #include "gridfire/engine/procedures/engine_procedures.h"
-#include "gridfire/engine/types/engine_types.h"
\ No newline at end of file
+#include "gridfire/engine/types/engine_types.h"
+#include "gridfire/engine/diagnostics/dynamic_engine_diagnostics.h"
\ No newline at end of file
diff --git a/src/include/gridfire/engine/procedures/construction.h b/src/include/gridfire/engine/procedures/construction.h
index e8297ba7..64963a44 100644
--- a/src/include/gridfire/engine/procedures/construction.h
+++ b/src/include/gridfire/engine/procedures/construction.h
@@ -24,23 +24,21 @@ namespace gridfire::engine {
     enum class NetworkConstructionFlags : uint32_t {
         NONE = 0,
 
-        STRONG = 1 << 0, // 1
+        REACLIB_STRONG = 1 << 0, // 1
 
-        BETA_MINUS = 1 << 1, // 2
-        BETA_PLUS = 1 << 2, // 4
-        ELECTRON_CAPTURE = 1 << 3, // 8
-        POSITRON_CAPTURE = 1 << 4, // 16
+        WRL_BETA_MINUS = 1 << 1, // 2
+        WRL_BETA_PLUS = 1 << 2, // 4
+        WRL_ELECTRON_CAPTURE = 1 << 3, // 8
+        WRL_POSITRON_CAPTURE = 1 << 4, // 16
 
         REACLIB_WEAK = 1 << 5,
 
-        WRL_WEAK = BETA_MINUS | BETA_PLUS | ELECTRON_CAPTURE | POSITRON_CAPTURE,
+        WRL_WEAK = WRL_BETA_MINUS | WRL_BETA_PLUS | WRL_ELECTRON_CAPTURE | WRL_POSITRON_CAPTURE,
 
-        REACLIB = STRONG | REACLIB_WEAK,
+        REACLIB = REACLIB_STRONG | REACLIB_WEAK,
 
         // Currently we default to just reaclib reactions but include both their strong and weak set
         DEFAULT = REACLIB,
-
-        ALL = STRONG | WRL_WEAK
     };
 
     /** @brief Helper function to convert NetworkConstructionFlags to their underlying integer type.
@@ -103,20 +101,20 @@ namespace gridfire::engine {
     inline std::string NetworkConstructionFlagsToString(NetworkConstructionFlags flags) {
         std::stringstream ss;
         constexpr std::array<NetworkConstructionFlags, 6> bases_flags_array = {
-            NetworkConstructionFlags::STRONG,
-            NetworkConstructionFlags::BETA_MINUS,
-            NetworkConstructionFlags::BETA_PLUS,
-            NetworkConstructionFlags::ELECTRON_CAPTURE,
-            NetworkConstructionFlags::POSITRON_CAPTURE,
+            NetworkConstructionFlags::REACLIB_STRONG,
+            NetworkConstructionFlags::WRL_BETA_MINUS,
+            NetworkConstructionFlags::WRL_BETA_PLUS,
+            NetworkConstructionFlags::WRL_ELECTRON_CAPTURE,
+            NetworkConstructionFlags::WRL_POSITRON_CAPTURE,
             NetworkConstructionFlags::REACLIB_WEAK
         };
 
         const std::unordered_map<NetworkConstructionFlags, std::string> bases_string_map = {
-            {NetworkConstructionFlags::STRONG, "Strong"},
-            {NetworkConstructionFlags::BETA_MINUS, "BetaMinus"},
-            {NetworkConstructionFlags::BETA_PLUS, "BetaPlus"},
-            {NetworkConstructionFlags::ELECTRON_CAPTURE, "ElectronCapture"},
-            {NetworkConstructionFlags::POSITRON_CAPTURE, "PositronCapture"},
+            {NetworkConstructionFlags::REACLIB_STRONG, "Strong"},
+            {NetworkConstructionFlags::WRL_BETA_MINUS, "BetaMinus"},
+            {NetworkConstructionFlags::WRL_BETA_PLUS, "BetaPlus"},
+            {NetworkConstructionFlags::WRL_ELECTRON_CAPTURE, "ElectronCapture"},
+            {NetworkConstructionFlags::WRL_POSITRON_CAPTURE, "PositronCapture"},
             {NetworkConstructionFlags::REACLIB_WEAK, "ReaclibWeak"}
         };
 
diff --git a/src/include/gridfire/engine/procedures/priming.h b/src/include/gridfire/engine/procedures/priming.h
index ed369668..65988270 100644
--- a/src/include/gridfire/engine/procedures/priming.h
+++ b/src/include/gridfire/engine/procedures/priming.h
@@ -31,54 +31,4 @@ namespace gridfire::engine {
          GraphEngine& engine,
          const std::optional<std::vector<reaction::ReactionType>>& ignoredReactionTypes
      );
-
-    /**
-     * @brief Computes the destruction rate constant for a specific species.
-     *
-     * Calculates the sum of molar reaction flows for all reactions where the species
-     * is a reactant (negative stoichiometry) after scaling its abundance to unity.
-     *
-     * @param engine Engine providing the current set of network reactions and flow calculations.
-     * @param species The atomic species whose destruction rate is computed.
-     * @param composition Current composition providing abundances for all species.
-     * @param T9 Temperature in units of 10^9 K.
-     * @param rho Density of the medium.
-     * @param reactionTypesToIgnore types of reactions to ignore during calculation.
-     * @pre Y.size() matches engine.getNetworkReactions().size() mapping species order.
-     * @post Returned rate constant is non-negative.
-     * @return Sum of absolute stoichiometry-weighted destruction flows for the species.
-     */
-     double calculateDestructionRateConstant(
-         const DynamicEngine& engine,
-         const fourdst::atomic::Species& species,
-         const fourdst::composition::Composition& composition,
-         double T9,
-         double rho,
-         const std::optional<std::vector<reaction::ReactionType>> &reactionTypesToIgnore
-     );
-
-    /**
-     * @brief Computes the creation rate for a specific species.
-     *
-     * Sums molar reaction flows for all reactions where the species
-     * appears as a product (positive stoichiometry).
-     *
-     * @param engine Engine providing the current set of network reactions and flow calculations.
-     * @param species The atomic species whose creation rate is computed.
-     * @param composition Composition object containing current abundances.
-     * @param T9 Temperature in units of 10^9 K.
-     * @param rho Density of the medium.
-     * @param reactionTypesToIgnore types of reactions to ignore during calculation.
-     * @pre Y.size() matches engine.getNetworkReactions().size() mapping species order.
-     * @post Returned creation rate is non-negative.
-     * @return Sum of stoichiometry-weighted creation flows for the species.
-     */
-     double calculateCreationRate(
-         const DynamicEngine& engine,
-         const fourdst::atomic::Species& species,
-         const fourdst::composition::Composition& composition,
-         double T9,
-         double rho,
-         const std::optional<std::vector<reaction::ReactionType>> &reactionTypesToIgnore
-     );
  }
\ No newline at end of file
diff --git a/src/include/gridfire/engine/types/engine_types.h b/src/include/gridfire/engine/types/engine_types.h
index 2baef958..f6d4b746 100644
--- a/src/include/gridfire/engine/types/engine_types.h
+++ b/src/include/gridfire/engine/types/engine_types.h
@@ -1,5 +1,7 @@
 #pragma once
 
+#include <string_view>
+
 namespace gridfire::engine {
     /**
      * @enum EngineTypes
diff --git a/src/include/gridfire/policy/policy_abstract.h b/src/include/gridfire/policy/policy_abstract.h
index c4f180ee..e89fb202 100644
--- a/src/include/gridfire/policy/policy_abstract.h
+++ b/src/include/gridfire/policy/policy_abstract.h
@@ -153,7 +153,7 @@ namespace gridfire::policy {
          * if (s != NetworkPolicyStatus::INITIALIZED_VERIFIED) { // handle error }
          * @endcode
          */
-        [[nodiscard]] virtual NetworkPolicyStatus getStatus() const = 0;
+        [[nodiscard]] virtual NetworkPolicyStatus get_status() const = 0;
 
         [[nodiscard]] virtual const std::vector<std::unique_ptr<engine::DynamicEngine>> &get_engine_stack() const = 0;
 
diff --git a/src/include/gridfire/policy/stellar_policy.h b/src/include/gridfire/policy/stellar_policy.h
index ddc0d499..89b2de39 100644
--- a/src/include/gridfire/policy/stellar_policy.h
+++ b/src/include/gridfire/policy/stellar_policy.h
@@ -141,7 +141,7 @@ namespace gridfire::policy {
          * @brief Gets the current status of the policy.
          * @return NetworkPolicyStatus The construction and verification status.
          */
-        [[nodiscard]] NetworkPolicyStatus getStatus() const override;
+        [[nodiscard]] NetworkPolicyStatus get_status() const override;
 
         [[nodiscard]] const std::vector<std::unique_ptr<engine::DynamicEngine>> &get_engine_stack() const override;
 
diff --git a/src/lib/engine/diagnostics/dynamic_engine_diagnostics.cpp b/src/lib/engine/diagnostics/dynamic_engine_diagnostics.cpp
index 3366247f..69719e09 100644
--- a/src/lib/engine/diagnostics/dynamic_engine_diagnostics.cpp
+++ b/src/lib/engine/diagnostics/dynamic_engine_diagnostics.cpp
@@ -15,7 +15,7 @@ namespace gridfire::engine::diagnostics {
         const double relTol,
         const double absTol,
         const size_t top_n,
-        bool json
+        const bool json
     ) {
         struct SpeciesError {
             std::string name;
diff --git a/src/lib/engine/procedures/construction.cpp b/src/lib/engine/procedures/construction.cpp
index 91778ea5..ec11fbe0 100644
--- a/src/lib/engine/procedures/construction.cpp
+++ b/src/lib/engine/procedures/construction.cpp
@@ -63,7 +63,7 @@ namespace {
                 parent_species.z() - 1
             );
             if (downProduct.has_value()) { // Only add the reaction if the Species map contains the product
-                if (has_flag(reactionTypes, gridfire::engine::NetworkConstructionFlags::BETA_PLUS)) {
+                if (has_flag(reactionTypes, gridfire::engine::NetworkConstructionFlags::WRL_BETA_PLUS)) {
                     weak_reaction_pool.add_reaction(
                         std::make_unique<gridfire::rates::weak::WeakReaction>(
                             parent_species,
@@ -72,7 +72,7 @@ namespace {
                         )
                     );
                 }
-                if (has_flag(reactionTypes, gridfire::engine::NetworkConstructionFlags::ELECTRON_CAPTURE)) {
+                if (has_flag(reactionTypes, gridfire::engine::NetworkConstructionFlags::WRL_ELECTRON_CAPTURE)) {
                     weak_reaction_pool.add_reaction(
                         std::make_unique<gridfire::rates::weak::WeakReaction>(
                             parent_species,
@@ -83,7 +83,7 @@ namespace {
                 }
             }
             if (upProduct.has_value()) { // Only add the reaction if the Species map contains the product
-                if (has_flag(reactionTypes, gridfire::engine::NetworkConstructionFlags::BETA_MINUS)) {
+                if (has_flag(reactionTypes, gridfire::engine::NetworkConstructionFlags::WRL_BETA_MINUS)) {
                     weak_reaction_pool.add_reaction(
                         std::make_unique<gridfire::rates::weak::WeakReaction>(
                             parent_species,
@@ -92,7 +92,7 @@ namespace {
                         )
                     );
                 }
-                if (has_flag(reactionTypes, gridfire::engine::NetworkConstructionFlags::POSITRON_CAPTURE)) {
+                if (has_flag(reactionTypes, gridfire::engine::NetworkConstructionFlags::WRL_POSITRON_CAPTURE)) {
                     weak_reaction_pool.add_reaction(
                         std::make_unique<gridfire::rates::weak::WeakReaction>(
                             parent_species,
@@ -111,7 +111,7 @@ namespace {
         const gridfire::engine::NetworkConstructionFlags reaction_types
     ) {
         gridfire::reaction::ReactionSet strong_reaction_pool;
-        if (has_flag(reaction_types, gridfire::engine::NetworkConstructionFlags::STRONG)) {
+        if (has_flag(reaction_types, gridfire::engine::NetworkConstructionFlags::REACLIB_STRONG)) {
             const auto& allReaclibReactions = gridfire::reaclib::get_all_reaclib_reactions();
             for (const auto& reaction : allReaclibReactions) {
                 const bool isWeakReaction = reaclib_reaction_is_weak(*reaction);
@@ -129,11 +129,11 @@ namespace {
     bool validate_unique_weak_set(gridfire::engine::NetworkConstructionFlags flag) {
         // This method ensures that weak reactions will only be fetched from either reaclib or the weak reaction library (WRL)
         // but not both
-        std::array<gridfire::engine::NetworkConstructionFlags, 4> WRL_Flags = {
-            gridfire::engine::NetworkConstructionFlags::BETA_PLUS,
-            gridfire::engine::NetworkConstructionFlags::ELECTRON_CAPTURE,
-            gridfire::engine::NetworkConstructionFlags::POSITRON_CAPTURE,
-            gridfire::engine::NetworkConstructionFlags::BETA_MINUS
+        const std::array<gridfire::engine::NetworkConstructionFlags, 4> WRL_Flags = {
+            gridfire::engine::NetworkConstructionFlags::WRL_BETA_PLUS,
+            gridfire::engine::NetworkConstructionFlags::WRL_ELECTRON_CAPTURE,
+            gridfire::engine::NetworkConstructionFlags::WRL_POSITRON_CAPTURE,
+            gridfire::engine::NetworkConstructionFlags::WRL_BETA_MINUS
         };
 
         if (!has_flag(flag, gridfire::engine::NetworkConstructionFlags::REACLIB_WEAK)) {
diff --git a/src/lib/policy/stellar_policy.cpp b/src/lib/policy/stellar_policy.cpp
index f17ce5d6..22dbcc32 100644
--- a/src/lib/policy/stellar_policy.cpp
+++ b/src/lib/policy/stellar_policy.cpp
@@ -92,7 +92,7 @@ namespace gridfire::policy {
         return std::make_unique<partition::CompositePartitionFunction>(partitionFunction);
     }
 
-    inline NetworkPolicyStatus MainSequencePolicy::getStatus() const {
+    inline NetworkPolicyStatus MainSequencePolicy::get_status() const {
         return m_status;
     }
 
diff --git a/src/python/bindings.cpp b/src/python/bindings.cpp
index 0d29b21b..b0d88baa 100644
--- a/src/python/bindings.cpp
+++ b/src/python/bindings.cpp
@@ -3,7 +3,6 @@
 
 #include "types/bindings.h"
 #include "partition/bindings.h"
-#include "expectations/bindings.h"
 #include "engine/bindings.h"
 #include "exceptions/bindings.h"
 #include "io/bindings.h"
@@ -11,8 +10,9 @@
 #include "screening/bindings.h"
 #include "solver/bindings.h"
 #include "utils/bindings.h"
+#include "policy/bindings.h"
 
-PYBIND11_MODULE(gridfire, m) {
+PYBIND11_MODULE(_gridfire, m) {
     m.doc() = "Python bindings for the fourdst utility modules which are a part of the 4D-STAR project.";
 
     pybind11::module::import("fourdst.constants");
@@ -26,9 +26,6 @@ PYBIND11_MODULE(gridfire, m) {
     auto partitionMod  = m.def_submodule("partition", "GridFire partition function bindings");
     register_partition_bindings(partitionMod);
 
-    auto expectationMod  = m.def_submodule("expectations", "GridFire expectations bindings");
-    register_expectation_bindings(expectationMod);
-
     auto reactionMod  = m.def_submodule("reaction", "GridFire reaction bindings");
     register_reaction_bindings(reactionMod);
 
@@ -47,6 +44,9 @@ PYBIND11_MODULE(gridfire, m) {
     auto solverMod  = m.def_submodule("solver", "GridFire numerical solver bindings");
     register_solver_bindings(solverMod);
 
+    auto policyMod  = m.def_submodule("policy", "GridFire network policy bindings");
+    register_policy_bindings(policyMod);
+
     auto utilsMod  = m.def_submodule("utils", "GridFire utility method bindings");
     register_utils_bindings(utilsMod);
 }
\ No newline at end of file
diff --git a/src/python/engine/bindings.cpp b/src/python/engine/bindings.cpp
index c89e8816..93675194 100644
--- a/src/python/engine/bindings.cpp
+++ b/src/python/engine/bindings.cpp
@@ -5,8 +5,8 @@
 #include "bindings.h"
 
 #include "gridfire/engine/engine.h"
-#include "gridfire/engine/diagnostics/dynamic_engine_diagnostics.h"
 #include "gridfire/exceptions/exceptions.h"
+#include "pybind11/numpy.h"
 #include "trampoline/py_engine.h"
 
 
@@ -15,21 +15,21 @@ namespace py = pybind11;
 
 namespace {
     template <typename T>
-    concept IsDynamicEngine = std::is_base_of_v<gridfire::DynamicEngine, T>;
+    concept IsDynamicEngine = std::is_base_of_v<gridfire::engine::DynamicEngine, T>;
 
     template <IsDynamicEngine T, IsDynamicEngine BaseT>
     void registerDynamicEngineDefs(py::class_<T, BaseT> pyClass) {
         pyClass.def(
             "calculateRHSAndEnergy",
             [](
-                const gridfire::DynamicEngine& self,
+                const gridfire::engine::DynamicEngine& self,
                 const fourdst::composition::Composition& comp,
                 const double T9,
                 const double rho
             ) {
                 auto result = self.calculateRHSAndEnergy(comp, T9, rho);
                 if (!result.has_value()) {
-                    throw gridfire::exceptions::StaleEngineError("Engine reports stale state, call update().");
+                    throw gridfire::exceptions::EngineError(std::format("calculateRHSAndEnergy returned a potentially recoverable error {}", gridfire::engine::EngineStatus_to_string(result.error())));
                 }
                 return result.value();
             },
@@ -39,21 +39,32 @@ namespace {
             "Calculate the right-hand side (dY/dt) and energy generation rate."
         )
         .def("calculateEpsDerivatives",
-            &gridfire::DynamicEngine::calculateEpsDerivatives,
+            &gridfire::engine::DynamicEngine::calculateEpsDerivatives,
             py::arg("comp"),
             py::arg("T9"),
             py::arg("rho"),
             "Calculate deps/dT and deps/drho"
         )
         .def("generateJacobianMatrix",
-            py::overload_cast<const fourdst::composition::Composition&, double, double>(&T::generateJacobianMatrix, py::const_),
+            [](const gridfire::engine::DynamicEngine& self,
+               const fourdst::composition::Composition& comp,
+               const double T9,
+               const double rho) -> gridfire::engine::NetworkJacobian {
+                return self.generateJacobianMatrix(comp, T9, rho);
+            },
             py::arg("comp"),
             py::arg("T9"),
             py::arg("rho"),
             "Generate the Jacobian matrix for the current state."
         )
         .def("generateJacobianMatrix",
-            py::overload_cast<const fourdst::composition::Composition&, double, double, const std::vector<fourdst::atomic::Species>&>(&T::generateJacobianMatrix, py::const_),
+            [](const gridfire::engine::DynamicEngine& self,
+               const fourdst::composition::Composition& comp,
+               const double T9,
+               const double rho,
+               const std::vector<fourdst::atomic::Species>& activeSpecies) -> gridfire::engine::NetworkJacobian {
+                return self.generateJacobianMatrix(comp, T9, rho, activeSpecies);
+            },
             py::arg("comp"),
             py::arg("T9"),
             py::arg("rho"),
@@ -61,7 +72,13 @@ namespace {
             "Generate the jacobian matrix only for the subset of the matrix representing the active species."
         )
         .def("generateJacobianMatrix",
-            py::overload_cast<const fourdst::composition::Composition&, double, double, const gridfire::SparsityPattern&>(&T::generateJacobianMatrix, py::const_),
+            [](const gridfire::engine::DynamicEngine& self,
+               const fourdst::composition::Composition& comp,
+               const double T9,
+               const double rho,
+               const gridfire::engine::SparsityPattern& sparsityPattern) -> gridfire::engine::NetworkJacobian {
+                return self.generateJacobianMatrix(comp, T9, rho, sparsityPattern);
+            },
             py::arg("comp"),
             py::arg("T9"),
             py::arg("rho"),
@@ -73,7 +90,7 @@ namespace {
         )
         .def("calculateMolarReactionFlow",
             [](
-                const gridfire::DynamicEngine& self,
+                const gridfire::engine::DynamicEngine& self,
                 const gridfire::reaction::Reaction& reaction,
                 const fourdst::composition::Composition& comp,
                 const double T9,
@@ -97,11 +114,6 @@ namespace {
             py::arg("reactions"),
             "Set the network reactions to a new set of reactions."
         )
-        .def("getJacobianMatrixEntry", &T::getJacobianMatrixEntry,
-            py::arg("rowSpecies"),
-            py::arg("colSpecies"),
-            "Get an entry from the previously generated Jacobian matrix."
-        )
         .def("getStoichiometryMatrixEntry", &T::getStoichiometryMatrixEntry,
             py::arg("species"),
             py::arg("reaction"),
@@ -109,14 +121,14 @@ namespace {
         )
         .def("getSpeciesTimescales",
             [](
-                const gridfire::DynamicEngine& self,
+                const gridfire::engine::DynamicEngine& self,
                 const fourdst::composition::Composition& comp,
                 const double T9,
                 const double rho
             ) -> std::unordered_map<fourdst::atomic::Species, double> {
                 const auto result = self.getSpeciesTimescales(comp, T9, rho);
                 if (!result.has_value()) {
-                    throw gridfire::exceptions::StaleEngineError("Engine reports stale state, call update().");
+                    throw gridfire::exceptions::EngineError(std::format("getSpeciesTimescales has returned a potentially recoverable error {}", gridfire::engine::EngineStatus_to_string(result.error())));
                 }
                 return result.value();
             },
@@ -127,14 +139,14 @@ namespace {
         )
         .def("getSpeciesDestructionTimescales",
             [](
-                const gridfire::DynamicEngine& self,
+                const gridfire::engine::DynamicEngine& self,
                 const fourdst::composition::Composition& comp,
                 const double T9,
                 const double rho
             ) -> std::unordered_map<fourdst::atomic::Species, double> {
                 const auto result = self.getSpeciesDestructionTimescales(comp, T9, rho);
                 if (!result.has_value()) {
-                    throw gridfire::exceptions::StaleEngineError("Engine reports stale state, call update().");
+                    throw gridfire::exceptions::EngineError(std::format("getSpeciesDestructionTimescales has returned a potentially recoverable error {}", gridfire::engine::EngineStatus_to_string(result.error())));
                 }
                 return result.value();
             },
@@ -179,7 +191,7 @@ namespace {
         .def("rebuild",
             &T::rebuild,
             py::arg("composition"),
-            py::arg("depth") = gridfire::NetworkBuildDepth::Full,
+            py::arg("depth") = gridfire::engine::NetworkBuildDepth::Full,
             "Rebuild the engine with a new composition and build depth."
         )
         .def("isStale",
@@ -190,7 +202,14 @@ namespace {
         .def("collectComposition",
             &T::collectComposition,
             py::arg("composition"),
+            py::arg("T9"),
+            py::arg("rho"),
             "Recursively collect composition from current engine and any sub engines if they exist."
+        )
+        .def("getSpeciesStatus",
+            &T::getSpeciesStatus,
+            py::arg("species"),
+            "Get the status of a species in the network."
         );
 
     }
@@ -206,11 +225,11 @@ void register_engine_bindings(py::module &m) {
 
 void register_base_engine_bindings(const pybind11::module &m) {
 
-    py::class_<gridfire::StepDerivatives<double>>(m, "StepDerivatives")
-        .def_readonly("dYdt", &gridfire::StepDerivatives<double>::dydt, "The right-hand side (dY/dt) of the ODE system.")
-        .def_readonly("energy", &gridfire::StepDerivatives<double>::nuclearEnergyGenerationRate, "The energy generation rate.");
+    py::class_<gridfire::engine::StepDerivatives<double>>(m, "StepDerivatives")
+        .def_readonly("dYdt", &gridfire::engine::StepDerivatives<double>::dydt, "The right-hand side (dY/dt) of the ODE system.")
+        .def_readonly("energy", &gridfire::engine::StepDerivatives<double>::nuclearEnergyGenerationRate, "The energy generation rate.");
 
-    py::class_<gridfire::SparsityPattern> py_sparsity_pattern(m, "SparsityPattern");
+    py::class_<gridfire::engine::SparsityPattern> py_sparsity_pattern(m, "SparsityPattern");
 
     abs_stype_register_engine_bindings(m);
     abs_stype_register_dynamic_engine_bindings(m);
@@ -218,148 +237,219 @@ void register_base_engine_bindings(const pybind11::module &m) {
 }
 
 void abs_stype_register_engine_bindings(const pybind11::module &m) {
-    py::class_<gridfire::Engine, PyEngine>(m, "Engine");
+    py::class_<gridfire::engine::Engine, PyEngine>(m, "Engine");
 }
 
 void abs_stype_register_dynamic_engine_bindings(const pybind11::module &m) {
-    const auto a = py::class_<gridfire::DynamicEngine, PyDynamicEngine>(m, "DynamicEngine");
+    const auto a = py::class_<gridfire::engine::DynamicEngine, PyDynamicEngine>(m, "DynamicEngine");
 }
 
 void register_engine_procedural_bindings(pybind11::module &m) {
-    auto procedures = m.def_submodule("procedures", "Procedural functions associated with engine module");
-    register_engine_construction_bindings(procedures);
-    register_engine_construction_bindings(procedures);
+    register_engine_construction_bindings(m);
+    register_engine_priming_bindings(m);
 }
 
 void register_engine_diagnostic_bindings(pybind11::module &m) {
     auto diagnostics = m.def_submodule("diagnostics", "A submodule for engine diagnostics");
     diagnostics.def("report_limiting_species",
-        &gridfire::diagnostics::report_limiting_species,
+        &gridfire::engine::diagnostics::report_limiting_species,
         py::arg("engine"),
         py::arg("Y_full"),
         py::arg("E_full"),
-        py::arg("dydt_full"),
         py::arg("relTol"),
         py::arg("absTol"),
-        py::arg("top_n") = 10
+        py::arg("top_n"),
+        py::arg("json")
     );
 
     diagnostics.def("inspect_species_balance",
-        &gridfire::diagnostics::inspect_species_balance,
+        &gridfire::engine::diagnostics::inspect_species_balance,
         py::arg("engine"),
         py::arg("species_name"),
         py::arg("comp"),
         py::arg("T9"),
-        py::arg("rho")
+        py::arg("rho"),
+        py::arg("json")
     );
 
     diagnostics.def("inspect_jacobian_stiffness",
-        &gridfire::diagnostics::inspect_jacobian_stiffness,
+        &gridfire::engine::diagnostics::inspect_jacobian_stiffness,
         py::arg("engine"),
         py::arg("comp"),
         py::arg("T9"),
-        py::arg("rho")
+        py::arg("rho"),
+        py::arg("json")
     );
 }
 
 void register_engine_construction_bindings(pybind11::module &m) {
-    m.def("build_nuclear_network", &gridfire::build_nuclear_network,
+    py::enum_<gridfire::engine::NetworkConstructionFlags>(m, "NetworkConstructionFlags")
+        .value("NONE", gridfire::engine::NetworkConstructionFlags::NONE, "No special construction flags.")
+        .value("REACLIB_STRONG", gridfire::engine::NetworkConstructionFlags::REACLIB_STRONG, "Include strong reactions from reaclib.")
+        .value("WRL_BETA_MINUS", gridfire::engine::NetworkConstructionFlags::WRL_BETA_MINUS, "Include beta-minus decay reactions from weak rate library.")
+        .value("WRL_BETA_PLUS", gridfire::engine::NetworkConstructionFlags::WRL_BETA_PLUS, "Include beta-plus decay reactions from weak rate library.")
+        .value("WRL_ELECTRON_CAPTURE", gridfire::engine::NetworkConstructionFlags::WRL_ELECTRON_CAPTURE, "Include electron capture reactions from weak rate library.")
+        .value("WRL_POSITRON_CAPTURE", gridfire::engine::NetworkConstructionFlags::WRL_POSITRON_CAPTURE, "Include positron capture reactions from weak rate library.")
+        .value("REACLIB_WEAK", gridfire::engine::NetworkConstructionFlags::REACLIB_WEAK, "Include weak reactions from reaclib.")
+        .value("WRL_WEAK", gridfire::engine::NetworkConstructionFlags::WRL_WEAK, "Include all weak reactions from weak rate library.")
+        .value("REACLIB", gridfire::engine::NetworkConstructionFlags::REACLIB, "Include all reactions from reaclib.")
+        .value("DEFAULT", gridfire::engine::NetworkConstructionFlags::DEFAULT, "Default construction flags (Reaclib strong and weak).")
+        .export_values()
+        .def("__repr__", [](const gridfire::engine::NetworkConstructionFlags& flags) {
+            return gridfire::engine::NetworkConstructionFlagsToString(flags);
+        });
+
+    m.def("build_nuclear_network", &gridfire::engine::build_nuclear_network,
         py::arg("composition"),
         py::arg("weakInterpolator"),
-        py::arg("maxLayers") = gridfire::NetworkBuildDepth::Full,
-        py::arg("reverse") = false,
+        py::arg("maxLayers") = gridfire::engine::NetworkBuildDepth::Full,
+        py::arg("ReactionTypes") = gridfire::engine::NetworkConstructionFlags::DEFAULT,
         "Build a nuclear network from a composition using all archived reaction data."
     );
 }
 
 void register_engine_priming_bindings(pybind11::module &m) {
-
-    m.def("calculateDestructionRateConstant",
-        &gridfire::calculateDestructionRateConstant,
+    m.def("primeNetwork",
+        &gridfire::engine::primeNetwork,
+        py::arg("netIn"),
         py::arg("engine"),
-        py::arg("species"),
-        py::arg("composition"),
-        py::arg("T9"),
-        py::arg("rho"),
-        py::arg("reactionTypesToIgnore")
-    );
-
-    m.def("calculateCreationRate",
-        &gridfire::calculateCreationRate,
-        py::arg("engine"),
-        py::arg("species"),
-        py::arg("composition"),
-        py::arg("T9"),
-        py::arg("rho"),
-        py::arg("reactionTypesToIgnore")
+        py::arg("ignoredReactionTypes") = std::nullopt,
+        "Prime a network with a short timescale ignition"
     );
 }
 
 void register_engine_type_bindings(pybind11::module &m) {
-    auto types = m.def_submodule("types", "Types associated with engine module");
-    register_engine_building_type_bindings(types);
-    register_engine_reporting_type_bindings(types);
+    register_engine_building_type_bindings(m);
+    register_engine_reporting_type_bindings(m);
+    register_engine_types_bindings(m);
+    register_jacobian_type_bindings(m);
 
 }
 
-void register_engine_building_type_bindings(pybind11::module &m) {
-    py::enum_<gridfire::NetworkBuildDepth>(m, "NetworkBuildDepth")
-        .value("Full", gridfire::NetworkBuildDepth::Full, "Full network build depth")
-        .value("Shallow", gridfire::NetworkBuildDepth::Shallow, "Shallow network build depth")
-        .value("SecondOrder", gridfire::NetworkBuildDepth::SecondOrder, "Second order network build depth")
-        .value("ThirdOrder", gridfire::NetworkBuildDepth::ThirdOrder, "Third order network build depth")
-        .value("FourthOrder", gridfire::NetworkBuildDepth::FourthOrder, "Fourth order network build depth")
-        .value("FifthOrder", gridfire::NetworkBuildDepth::FifthOrder, "Fifth order network build depth")
+void register_engine_building_type_bindings(const pybind11::module &m) {
+    py::enum_<gridfire::engine::NetworkBuildDepth>(m, "NetworkBuildDepth")
+        .value("Full", gridfire::engine::NetworkBuildDepth::Full, "Full network build depth")
+        .value("Shallow", gridfire::engine::NetworkBuildDepth::Shallow, "Shallow network build depth")
+        .value("SecondOrder", gridfire::engine::NetworkBuildDepth::SecondOrder, "Second order network build depth")
+        .value("ThirdOrder", gridfire::engine::NetworkBuildDepth::ThirdOrder, "Third order network build depth")
+        .value("FourthOrder", gridfire::engine::NetworkBuildDepth::FourthOrder, "Fourth order network build depth")
+        .value("FifthOrder", gridfire::engine::NetworkBuildDepth::FifthOrder, "Fifth order network build depth")
         .export_values();
 
-    py::class_<gridfire::BuildDepthType> py_build_depth_type(m, "BuildDepthType");
+    py::class_<gridfire::engine::BuildDepthType> py_build_depth_type(m, "BuildDepthType");
 }
 
-void register_engine_reporting_type_bindings(pybind11::module &m) {
-    py::enum_<gridfire::PrimingReportStatus>(m, "PrimingReportStatus")
-        .value("FULL_SUCCESS", gridfire::PrimingReportStatus::FULL_SUCCESS, "Priming was full successful.")
-        .value("NO_SPECIES_TO_PRIME", gridfire::PrimingReportStatus::NO_SPECIES_TO_PRIME, "No species to prime.")
-        .value("MAX_ITERATIONS_REACHED", gridfire::PrimingReportStatus::MAX_ITERATIONS_REACHED, "Maximum iterations reached during priming.")
-        .value("FAILED_TO_FINALIZE_COMPOSITION", gridfire::PrimingReportStatus::FAILED_TO_FINALIZE_COMPOSITION, "Failed to finalize the composition after priming.")
-        .value("FAILED_TO_FIND_CREATION_CHANNEL", gridfire::PrimingReportStatus::FAILED_TO_FIND_CREATION_CHANNEL, "Failed to find a creation channel for the priming species.")
-        .value("FAILED_TO_FIND_PRIMING_REACTIONS", gridfire::PrimingReportStatus::FAILED_TO_FIND_PRIMING_REACTIONS, "Failed to find priming reactions for the species.")
-        .value("BASE_NETWORK_TOO_SHALLOW", gridfire::PrimingReportStatus::BASE_NETWORK_TOO_SHALLOW, "The base network is too shallow for priming.")
+void register_engine_reporting_type_bindings(const pybind11::module &m) {
+    py::enum_<gridfire::engine::PrimingReportStatus>(m, "PrimingReportStatus")
+        .value("FULL_SUCCESS", gridfire::engine::PrimingReportStatus::SUCCESS, "Priming was full successful.")
+        .value("NO_SPECIES_TO_PRIME", gridfire::engine::PrimingReportStatus::SOLVER_FAILURE, "Solver Failed to converge during priming.")
+        .value("MAX_ITERATIONS_REACHED", gridfire::engine::PrimingReportStatus::ALREADY_PRIMED, "Engine has already been primed.")
         .export_values()
-        .def("__repr__", [](const gridfire::PrimingReportStatus& status) {
+        .def("__repr__", [](const gridfire::engine::PrimingReportStatus& status) {
                 std::stringstream ss;
-                ss << gridfire::PrimingReportStatusStrings.at(status) << "\n";
+                ss << gridfire::engine::PrimingReportStatusStrings.at(status) << "\n";
                 return ss.str();
             },
         "String representation of the PrimingReport."
         );
 
-    py::class_<gridfire::PrimingReport>(m, "PrimingReport")
-        .def_readonly("success", &gridfire::PrimingReport::success, "Indicates if the priming was successful.")
-        .def_readonly("massFractionChanges", &gridfire::PrimingReport::massFractionChanges, "Map of species to their mass fraction changes after priming.")
-        .def_readonly("primedComposition", &gridfire::PrimingReport::primedComposition, "The composition after priming.")
-        .def_readonly("status", &gridfire::PrimingReport::status, "Status message from the priming process.")
-        .def("__repr__", [](const gridfire::PrimingReport& report) {
+    py::class_<gridfire::engine::PrimingReport>(m, "PrimingReport")
+        .def_readonly("success", &gridfire::engine::PrimingReport::success, "Indicates if the priming was successful.")
+        .def_readonly("primedComposition", &gridfire::engine::PrimingReport::primedComposition, "The composition after priming.")
+        .def_readonly("status", &gridfire::engine::PrimingReport::status, "Status message from the priming process.")
+        .def("__repr__", [](const gridfire::engine::PrimingReport& report) {
                 std::stringstream ss;
                 ss << report;
                 return ss.str();
             }
         );
+
+    py::enum_<gridfire::engine::SpeciesStatus>(m, "SpeciesStatus")
+        .value("ACTIVE", gridfire::engine::SpeciesStatus::ACTIVE, "Species is active in the network.")
+        .value("EQUILIBRIUM", gridfire::engine::SpeciesStatus::EQUILIBRIUM, "Species is in equilibrium.")
+        .value("INACTIVE_FLOW", gridfire::engine::SpeciesStatus::INACTIVE_FLOW, "Species is inactive due to flow.")
+        .value("NOT_PRESENT", gridfire::engine::SpeciesStatus::NOT_PRESENT, "Species is not present in the network.")
+        .export_values()
+        .def("__repr__", [](const gridfire::engine::SpeciesStatus& status) {
+            return gridfire::engine::SpeciesStatus_to_string(status);
+        });
 }
 
+void register_engine_types_bindings(const pybind11::module &m) {
+    py::enum_<gridfire::engine::EngineTypes>(m, "EngineTypes")
+        .value("GRAPH_ENGINE", gridfire::engine::EngineTypes::GRAPH_ENGINE, "The standard graph-based engine.")
+        .value("ADAPTIVE_ENGINE_VIEW", gridfire::engine::EngineTypes::ADAPTIVE_ENGINE_VIEW, "An engine that adapts based on certain criteria.")
+        .value("MULTISCALE_PARTITIONING_ENGINE_VIEW", gridfire::engine::EngineTypes::MULTISCALE_PARTITIONING_ENGINE_VIEW, "An engine that partitions the system at multiple scales.")
+        .value("PRIMING_ENGINE_VIEW", gridfire::engine::EngineTypes::PRIMING_ENGINE_VIEW, "An engine that uses a priming strategy for simulations.")
+        .value("DEFINED_ENGINE_VIEW", gridfire::engine::EngineTypes::DEFINED_ENGINE_VIEW, "An engine defined by user specifications.")
+        .value("FILE_DEFINED_ENGINE_VIEW", gridfire::engine::EngineTypes::FILE_DEFINED_ENGINE_VIEW, "An engine defined through external files.")
+        .export_values()
+        .def("__repr__", [](const gridfire::engine::EngineTypes& type) {
+                return std::string(gridfire::engine::engine_type_to_string(type));
+            },
+        "String representation of the EngineTypes."
+        );
+}
+
+void register_jacobian_type_bindings(pybind11::module &m) {
+    py::class_<gridfire::engine::NetworkJacobian>(m, "NetworkJacobian")
+        .def("rank", &gridfire::engine::NetworkJacobian::rank, "Get the rank of the Jacobian matrix.")
+        .def("nnz", &gridfire::engine::NetworkJacobian::nnz, "Get the number of non-zero entries in the Jacobian matrix.")
+        .def("singular", &gridfire::engine::NetworkJacobian::singular, "Check if the Jacobian matrix is singular.")
+        .def("infs", &gridfire::engine::NetworkJacobian::infs, "Get all infinite entries in the Jacobian matrix.")
+        .def("nans", &gridfire::engine::NetworkJacobian::nans, "Get all NaN entries in the Jacobian matrix.")
+        .def("data", &gridfire::engine::NetworkJacobian::data, "Get the underlying sparse matrix data.")
+        .def("mapping", &gridfire::engine::NetworkJacobian::mapping, "Get the species-to-index mapping.")
+        .def("shape", &gridfire::engine::NetworkJacobian::shape, "Get the shape of the Jacobian matrix as (rows, columns).")
+        .def("to_csv", &gridfire::engine::NetworkJacobian::to_csv, py::arg("filename"), "Export the Jacobian matrix to a CSV file.")
+        .def("__getitem__", [](const gridfire::engine::NetworkJacobian &self, const std::pair<fourdst::atomic::Species, fourdst::atomic::Species>& speciesPair) {
+            return self(speciesPair.first, speciesPair.second);
+        }, py::arg("key"), "Get an entry from the Jacobian matrix using species identifiers.")
+        .def("__getitem__", [](const gridfire::engine::NetworkJacobian &self, const std::pair<size_t, size_t>& indexPair) {
+            return self(indexPair.first, indexPair.second);
+        }, py::arg("key"), "Get an entry from the Jacobian matrix using indices.")
+        .def("__setitem__", [](gridfire::engine::NetworkJacobian &self, const std::pair<fourdst::atomic::Species, fourdst::atomic::Species>& speciesPair, double value) {
+            self.set(speciesPair.first, speciesPair.second, value);
+        }, py::arg("key"), py::arg("value"), "Set an entry in the Jacobian matrix using species identifiers.")
+        .def("__setitem__", [](gridfire::engine::NetworkJacobian &self, const std::pair<size_t, size_t>& indexPair, double value) {
+            self.set(indexPair.first, indexPair.second, value);
+        }, py::arg("key"), py::arg("value"), "Set an entry in the Jacobian matrix using indices.")
+        .def("to_numpy", [](const gridfire::engine::NetworkJacobian &self) {
+            auto denseMatrix = Eigen::MatrixXd(self.data());
+            return pybind11::array_t<double>(
+                {std::get<0>(self.shape()), std::get<1>(self.shape())},
+                {sizeof(double) * std::get<1>(self.shape()), sizeof(double)},
+                denseMatrix.data()
+            );
+        }, "Convert the Jacobian matrix to a NumPy array.");
+
+    m.def(
+        "regularize_jacobian",
+        [](const gridfire::engine::NetworkJacobian& jac, const fourdst::composition::Composition& comp) {
+            return regularize_jacobian(jac, comp, std::nullopt);
+        },
+        py::arg("jacobian"),
+        py::arg("composition"),
+        "regularize_jacobian"
+    );
+}
+
+
+
 void con_stype_register_graph_engine_bindings(const pybind11::module &m) {
 
-    auto py_graph_engine_bindings = py::class_<gridfire::GraphEngine, gridfire::DynamicEngine>(m, "GraphEngine");
+    auto py_graph_engine_bindings = py::class_<gridfire::engine::GraphEngine, gridfire::engine::DynamicEngine>(m, "GraphEngine");
 
     // Register the Graph Engine Specific Bindings
-    py_graph_engine_bindings.def(py::init<const fourdst::composition::Composition &, const gridfire::BuildDepthType>(),
+    py_graph_engine_bindings.def(py::init<const fourdst::composition::Composition &, const gridfire::engine::BuildDepthType>(),
         py::arg("composition"),
-        py::arg("depth") = gridfire::NetworkBuildDepth::Full,
+        py::arg("depth") = gridfire::engine::NetworkBuildDepth::Full,
         "Initialize GraphEngine with a composition and build depth."
     );
-    py_graph_engine_bindings.def(py::init<const fourdst::composition::Composition &,const gridfire::partition::PartitionFunction &, const gridfire::BuildDepthType>(),
+    py_graph_engine_bindings.def(py::init<const fourdst::composition::Composition &,const gridfire::partition::PartitionFunction &, const gridfire::engine::BuildDepthType>(),
         py::arg("composition"),
         py::arg("partitionFunction"),
-        py::arg("depth") = gridfire::NetworkBuildDepth::Full,
+        py::arg("depth") = gridfire::engine::NetworkBuildDepth::Full,
         "Initialize GraphEngine with a composition, partition function and build depth."
     );
     py_graph_engine_bindings.def(py::init<const gridfire::reaction::ReactionSet &>(),
@@ -367,54 +457,66 @@ void con_stype_register_graph_engine_bindings(const pybind11::module &m) {
         "Initialize GraphEngine with a set of reactions."
     );
     py_graph_engine_bindings.def_static("getNetReactionStoichiometry",
-        &gridfire::GraphEngine::getNetReactionStoichiometry,
+        &gridfire::engine::GraphEngine::getNetReactionStoichiometry,
         py::arg("reaction"),
         "Get the net stoichiometry for a given reaction."
     );
     py_graph_engine_bindings.def("getSpeciesTimescales",
-        py::overload_cast<const fourdst::composition::Composition&, double, double, const gridfire::reaction::ReactionSet&>(&gridfire::GraphEngine::getSpeciesTimescales, py::const_),
+        [](const gridfire::engine::GraphEngine& self,
+           const fourdst::composition::Composition& composition,
+           const double T9,
+           const double rho,
+           const gridfire::reaction::ReactionSet& activeReactions) {
+            return self.getSpeciesTimescales(composition, T9, rho, activeReactions);
+        },
         py::arg("composition"),
         py::arg("T9"),
         py::arg("rho"),
         py::arg("activeReactions")
     );
     py_graph_engine_bindings.def("getSpeciesDestructionTimescales",
-        py::overload_cast<const fourdst::composition::Composition&, double, double, const gridfire::reaction::ReactionSet&>(&gridfire::GraphEngine::getSpeciesDestructionTimescales, py::const_),
+        [](const gridfire::engine::GraphEngine& self,
+           const fourdst::composition::Composition& composition,
+           const double T9,
+           const double rho,
+           const gridfire::reaction::ReactionSet& activeReactions) {
+            return self.getSpeciesDestructionTimescales(composition, T9, rho, activeReactions);
+        },
         py::arg("composition"),
         py::arg("T9"),
         py::arg("rho"),
         py::arg("activeReactions")
     );
     py_graph_engine_bindings.def("involvesSpecies",
-        &gridfire::GraphEngine::involvesSpecies,
+        &gridfire::engine::GraphEngine::involvesSpecies,
         py::arg("species"),
         "Check if a given species is involved in the network."
     );
     py_graph_engine_bindings.def("exportToDot",
-        &gridfire::GraphEngine::exportToDot,
+        &gridfire::engine::GraphEngine::exportToDot,
         py::arg("filename"),
         "Export the network to a DOT file for visualization."
     );
     py_graph_engine_bindings.def("exportToCSV",
-        &gridfire::GraphEngine::exportToCSV,
+        &gridfire::engine::GraphEngine::exportToCSV,
         py::arg("filename"),
         "Export the network to a CSV file for analysis."
     );
     py_graph_engine_bindings.def("setPrecomputation",
-        &gridfire::GraphEngine::setPrecomputation,
+        &gridfire::engine::GraphEngine::setPrecomputation,
         py::arg("precompute"),
         "Enable or disable precomputation for the engine."
     );
     py_graph_engine_bindings.def("isPrecomputationEnabled",
-        &gridfire::GraphEngine::isPrecomputationEnabled,
+        &gridfire::engine::GraphEngine::isPrecomputationEnabled,
         "Check if precomputation is enabled for the engine."
     );
     py_graph_engine_bindings.def("getPartitionFunction",
-        &gridfire::GraphEngine::getPartitionFunction,
+        &gridfire::engine::GraphEngine::getPartitionFunction,
         "Get the partition function used by the engine."
     );
     py_graph_engine_bindings.def("calculateReverseRate",
-        &gridfire::GraphEngine::calculateReverseRate,
+        &gridfire::engine::GraphEngine::calculateReverseRate,
         py::arg("reaction"),
         py::arg("T9"),
         py::arg("rho"),
@@ -422,7 +524,7 @@ void con_stype_register_graph_engine_bindings(const pybind11::module &m) {
         "Calculate the reverse rate for a given reaction at a specific temperature, density, and composition."
     );
     py_graph_engine_bindings.def("calculateReverseRateTwoBody",
-        &gridfire::GraphEngine::calculateReverseRateTwoBody,
+        &gridfire::engine::GraphEngine::calculateReverseRateTwoBody,
         py::arg("reaction"),
         py::arg("T9"),
         py::arg("forwardRate"),
@@ -430,7 +532,7 @@ void con_stype_register_graph_engine_bindings(const pybind11::module &m) {
         "Calculate the reverse rate for a two-body reaction at a specific temperature."
     );
     py_graph_engine_bindings.def("calculateReverseRateTwoBodyDerivative",
-        &gridfire::GraphEngine::calculateReverseRateTwoBodyDerivative,
+        &gridfire::engine::GraphEngine::calculateReverseRateTwoBodyDerivative,
         py::arg("reaction"),
         py::arg("T9"),
         py::arg("rho"),
@@ -439,126 +541,96 @@ void con_stype_register_graph_engine_bindings(const pybind11::module &m) {
         "Calculate the derivative of the reverse rate for a two-body reaction at a specific temperature."
     );
     py_graph_engine_bindings.def("isUsingReverseReactions",
-        &gridfire::GraphEngine::isUsingReverseReactions,
+        &gridfire::engine::GraphEngine::isUsingReverseReactions,
         "Check if the engine is using reverse reactions."
     );
     py_graph_engine_bindings.def("setUseReverseReactions",
-        &gridfire::GraphEngine::setUseReverseReactions,
+        &gridfire::engine::GraphEngine::setUseReverseReactions,
         py::arg("useReverse"),
         "Enable or disable the use of reverse reactions in the engine."
     );
 
     // Register the general dynamic engine bindings
-    registerDynamicEngineDefs<gridfire::GraphEngine, gridfire::DynamicEngine>(py_graph_engine_bindings);
+    registerDynamicEngineDefs<gridfire::engine::GraphEngine, gridfire::engine::DynamicEngine>(py_graph_engine_bindings);
 }
 
 void register_engine_view_bindings(const pybind11::module &m) {
-    auto py_defined_engine_view_bindings = py::class_<gridfire::DefinedEngineView, gridfire::DynamicEngine>(m, "DefinedEngineView");
+    auto py_defined_engine_view_bindings = py::class_<gridfire::engine::DefinedEngineView, gridfire::engine::DynamicEngine>(m, "DefinedEngineView");
 
-    py_defined_engine_view_bindings.def(py::init<std::vector<std::string>, gridfire::GraphEngine&>(),
+    py_defined_engine_view_bindings.def(py::init<std::vector<std::string>, gridfire::engine::GraphEngine&>(),
         py::arg("peNames"),
         py::arg("baseEngine"),
         "Construct a defined engine view with a list of tracked reactions and a base engine."
     );
-    py_defined_engine_view_bindings.def("getBaseEngine", &gridfire::DefinedEngineView::getBaseEngine,
+    py_defined_engine_view_bindings.def("getBaseEngine", &gridfire::engine::DefinedEngineView::getBaseEngine,
         "Get the base engine associated with this defined engine view.");
 
-    registerDynamicEngineDefs<gridfire::DefinedEngineView, gridfire::DynamicEngine>(py_defined_engine_view_bindings);
+    registerDynamicEngineDefs<gridfire::engine::DefinedEngineView, gridfire::engine::DynamicEngine>(py_defined_engine_view_bindings);
 
-    auto py_file_defined_engine_view_bindings = py::class_<gridfire::FileDefinedEngineView, gridfire::DefinedEngineView>(m, "FileDefinedEngineView");
+    auto py_file_defined_engine_view_bindings = py::class_<gridfire::engine::FileDefinedEngineView, gridfire::engine::DefinedEngineView>(m, "FileDefinedEngineView");
     py_file_defined_engine_view_bindings.def(
-    py::init<gridfire::GraphEngine&, const std::string&, const gridfire::io::NetworkFileParser&>(),
+    py::init<gridfire::engine::GraphEngine&, const std::string&, const gridfire::io::NetworkFileParser&>(),
         py::arg("baseEngine"),
         py::arg("fileName"),
         py::arg("parser"),
         "Construct a defined engine view from a file and a base engine."
     );
-    py_file_defined_engine_view_bindings.def("getNetworkFile", &gridfire::FileDefinedEngineView::getNetworkFile,
+    py_file_defined_engine_view_bindings.def("getNetworkFile", &gridfire::engine::FileDefinedEngineView::getNetworkFile,
         "Get the network file associated with this defined engine view."
     );
-    py_file_defined_engine_view_bindings.def("getParser", &gridfire::FileDefinedEngineView::getParser,
+    py_file_defined_engine_view_bindings.def("getParser", &gridfire::engine::FileDefinedEngineView::getParser,
         "Get the parser used for this defined engine view."
     );
-    py_file_defined_engine_view_bindings.def("getBaseEngine", &gridfire::FileDefinedEngineView::getBaseEngine,
+    py_file_defined_engine_view_bindings.def("getBaseEngine", &gridfire::engine::FileDefinedEngineView::getBaseEngine,
         "Get the base engine associated with this file defined engine view.");
 
-    registerDynamicEngineDefs<gridfire::FileDefinedEngineView, gridfire::DefinedEngineView>(py_file_defined_engine_view_bindings);
+    registerDynamicEngineDefs<gridfire::engine::FileDefinedEngineView, gridfire::engine::DefinedEngineView>(py_file_defined_engine_view_bindings);
 
-    auto py_priming_engine_view_bindings = py::class_<gridfire::NetworkPrimingEngineView, gridfire::DefinedEngineView>(m, "NetworkPrimingEngineView");
-    py_priming_engine_view_bindings.def(py::init<const std::string&, gridfire::GraphEngine&>(),
+    auto py_priming_engine_view_bindings = py::class_<gridfire::engine::NetworkPrimingEngineView, gridfire::engine::DefinedEngineView>(m, "NetworkPrimingEngineView");
+    py_priming_engine_view_bindings.def(py::init<const std::string&, gridfire::engine::GraphEngine&>(),
         py::arg("primingSymbol"),
         py::arg("baseEngine"),
         "Construct a priming engine view with a priming symbol and a base engine.");
-    py_priming_engine_view_bindings.def(py::init<const fourdst::atomic::Species&, gridfire::GraphEngine&>(),
+    py_priming_engine_view_bindings.def(py::init<const fourdst::atomic::Species&, gridfire::engine::GraphEngine&>(),
         py::arg("primingSpecies"),
         py::arg("baseEngine"),
         "Construct a priming engine view with a priming species and a base engine.");
-    py_priming_engine_view_bindings.def("getBaseEngine", &gridfire::NetworkPrimingEngineView::getBaseEngine,
+    py_priming_engine_view_bindings.def("getBaseEngine", &gridfire::engine::NetworkPrimingEngineView::getBaseEngine,
         "Get the base engine associated with this priming engine view.");
 
-    registerDynamicEngineDefs<gridfire::NetworkPrimingEngineView, gridfire::DefinedEngineView>(py_priming_engine_view_bindings);
+    registerDynamicEngineDefs<gridfire::engine::NetworkPrimingEngineView, gridfire::engine::DefinedEngineView>(py_priming_engine_view_bindings);
 
-    auto py_adaptive_engine_view_bindings = py::class_<gridfire::AdaptiveEngineView, gridfire::DynamicEngine>(m, "AdaptiveEngineView");
-    py_adaptive_engine_view_bindings.def(py::init<gridfire::DynamicEngine&>(),
+    auto py_adaptive_engine_view_bindings = py::class_<gridfire::engine::AdaptiveEngineView, gridfire::engine::DynamicEngine>(m, "AdaptiveEngineView");
+    py_adaptive_engine_view_bindings.def(py::init<gridfire::engine::DynamicEngine&>(),
         py::arg("baseEngine"),
         "Construct an adaptive engine view with a base engine.");
     py_adaptive_engine_view_bindings.def("getBaseEngine",
-        &gridfire::AdaptiveEngineView::getBaseEngine,
+        &gridfire::engine::AdaptiveEngineView::getBaseEngine,
         "Get the base engine associated with this adaptive engine view."
     );
 
-    registerDynamicEngineDefs<gridfire::AdaptiveEngineView, gridfire::DynamicEngine>(py_adaptive_engine_view_bindings);
+    registerDynamicEngineDefs<gridfire::engine::AdaptiveEngineView, gridfire::engine::DynamicEngine>(py_adaptive_engine_view_bindings);
 
-    auto py_qse_cache_config = py::class_<gridfire::QSECacheConfig>(m, "QSECacheConfig");
-    auto py_qse_cache_key = py::class_<gridfire::QSECacheKey>(m, "QSECacheKey");
-
-    py_qse_cache_key.def(py::init<double, double, const std::vector<double>&>(),
-        py::arg("T9"),
-        py::arg("rho"),
-        py::arg("Y")
-    );
-
-    py_qse_cache_key.def("hash", &gridfire::QSECacheKey::hash,
-        "Get the pre-computed hash value of the key");
-
-    py_qse_cache_key.def_static("bin", &gridfire::QSECacheKey::bin,
-        py::arg("value"),
-        py::arg("tol"),
-        "bin a value based on a tolerance");
-    py_qse_cache_key.def("__eq__", &gridfire::QSECacheKey::operator==,
-        py::arg("other"),
-        "Check if two QSECacheKeys are equal");
-
-    auto py_multiscale_engine_view_bindings = py::class_<gridfire::MultiscalePartitioningEngineView, gridfire::DynamicEngine>(m, "MultiscalePartitioningEngineView");
-    py_multiscale_engine_view_bindings.def(py::init<gridfire::GraphEngine&>(),
+    auto py_multiscale_engine_view_bindings = py::class_<gridfire::engine::MultiscalePartitioningEngineView, gridfire::engine::DynamicEngine>(m, "MultiscalePartitioningEngineView");
+    py_multiscale_engine_view_bindings.def(py::init<gridfire::engine::GraphEngine&>(),
         py::arg("baseEngine"),
         "Construct a multiscale partitioning engine view with a base engine."
     );
     py_multiscale_engine_view_bindings.def("getBaseEngine",
-        &gridfire::MultiscalePartitioningEngineView::getBaseEngine,
+        &gridfire::engine::MultiscalePartitioningEngineView::getBaseEngine,
         "Get the base engine associated with this multiscale partitioning engine view."
     );
-    py_multiscale_engine_view_bindings.def("analyzeTimescalePoolConnectivity",
-        &gridfire::MultiscalePartitioningEngineView::analyzeTimescalePoolConnectivity,
-        py::arg("timescale_pools"),
-        py::arg("comp"),
-        py::arg("T9"),
-        py::arg("rho"),
-        "Analyze the connectivity of timescale pools in the network."
-    );
     py_multiscale_engine_view_bindings.def("partitionNetwork",
-        py::overload_cast<const fourdst::composition::Composition&, double, double>(&gridfire::MultiscalePartitioningEngineView::partitionNetwork),
-        py::arg("comp"),
-        py::arg("T9"),
-        py::arg("rho"),
+        &gridfire::engine::MultiscalePartitioningEngineView::partitionNetwork,
+        py::arg("netIn"),
         "Partition the network based on species timescales and connectivity.");
     py_multiscale_engine_view_bindings.def("partitionNetwork",
-        py::overload_cast<const gridfire::NetIn&>(&gridfire::MultiscalePartitioningEngineView::partitionNetwork),
+        py::overload_cast<const gridfire::NetIn&>(&gridfire::engine::MultiscalePartitioningEngineView::partitionNetwork),
         py::arg("netIn"),
         "Partition the network based on a NetIn object."
     );
     py_multiscale_engine_view_bindings.def("exportToDot",
-    &gridfire::MultiscalePartitioningEngineView::exportToDot,
+    &gridfire::engine::MultiscalePartitioningEngineView::exportToDot,
         py::arg("filename"),
         py::arg("comp"),
         py::arg("T9"),
@@ -566,26 +638,37 @@ void register_engine_view_bindings(const pybind11::module &m) {
         "Export the network to a DOT file for visualization."
     );
     py_multiscale_engine_view_bindings.def("getFastSpecies",
-        &gridfire::MultiscalePartitioningEngineView::getFastSpecies,
+        &gridfire::engine::MultiscalePartitioningEngineView::getFastSpecies,
         "Get the list of fast species in the network."
     );
     py_multiscale_engine_view_bindings.def("getDynamicSpecies",
-        &gridfire::MultiscalePartitioningEngineView::getDynamicSpecies,
+        &gridfire::engine::MultiscalePartitioningEngineView::getDynamicSpecies,
         "Get the list of dynamic species in the network."
     );
-    py_multiscale_engine_view_bindings.def("equilibrateNetwork",
-        py::overload_cast<const fourdst::composition::Composition&, double, double>(&gridfire::MultiscalePartitioningEngineView::equilibrateNetwork),
+    py_multiscale_engine_view_bindings.def("involvesSpecies",
+        &gridfire::engine::MultiscalePartitioningEngineView::involvesSpecies,
+        py::arg("species"),
+        "Check if a given species is involved in the network (in either the algebraic or dynamic set)."
+    );
+    py_multiscale_engine_view_bindings.def("involvesSpeciesInQSE",
+        &gridfire::engine::MultiscalePartitioningEngineView::involvesSpeciesInQSE,
+        py::arg("species"),
+        "Check if a given species is involved in the network's algebraic set."
+    );
+    py_multiscale_engine_view_bindings.def("involvesSpeciesInDynamic",
+        &gridfire::engine::MultiscalePartitioningEngineView::involvesSpeciesInDynamic,
+        py::arg("species"),
+        "Check if a given species is involved in the network's dynamic set."
+    );
+    py_multiscale_engine_view_bindings.def("getNormalizedEquilibratedComposition",
+        &gridfire::engine::MultiscalePartitioningEngineView::getNormalizedEquilibratedComposition,
         py::arg("comp"),
         py::arg("T9"),
         py::arg("rho"),
-        "Equilibrate the network based on species abundances and conditions.");
-    py_multiscale_engine_view_bindings.def("equilibrateNetwork",
-        py::overload_cast<const gridfire::NetIn&>(&gridfire::MultiscalePartitioningEngineView::equilibrateNetwork),
-        py::arg("netIn"),
-        "Equilibrate the network based on a NetIn object."
+        "Get the normalized equilibrated composition for the algebraic species."
     );
 
-    registerDynamicEngineDefs<gridfire::MultiscalePartitioningEngineView, gridfire::DynamicEngine>(
+    registerDynamicEngineDefs<gridfire::engine::MultiscalePartitioningEngineView, gridfire::engine::DynamicEngine>(
         py_multiscale_engine_view_bindings
     );
 
diff --git a/src/python/engine/bindings.h b/src/python/engine/bindings.h
index 86e0aac6..27b2454d 100644
--- a/src/python/engine/bindings.h
+++ b/src/python/engine/bindings.h
@@ -20,7 +20,9 @@ void register_engine_construction_bindings(pybind11::module &m);
 void register_engine_priming_bindings(pybind11::module &m);
 
 void register_engine_type_bindings(pybind11::module &m);
-void register_engine_building_type_bindings(pybind11::module &m);
-void register_engine_reporting_type_bindings(pybind11::module &m);
+void register_engine_building_type_bindings(const pybind11::module &m);
+void register_engine_reporting_type_bindings(const pybind11::module &m);
+void register_engine_types_bindings(const pybind11::module &m);
+void register_jacobian_type_bindings(pybind11::module &m);
 
 
diff --git a/src/python/engine/trampoline/py_engine.cpp b/src/python/engine/trampoline/py_engine.cpp
index abf25666..76da9ea4 100644
--- a/src/python/engine/trampoline/py_engine.cpp
+++ b/src/python/engine/trampoline/py_engine.cpp
@@ -33,10 +33,14 @@ const std::vector<fourdst::atomic::Species>& PyEngine::getNetworkSpecies() const
     py::pybind11_fail("Tried to call pure virtual function \"DynamicEngine::getNetworkSpecies\"");
 }
 
-std::expected<gridfire::StepDerivatives<double>, gridfire::expectations::StaleEngineError> PyEngine::calculateRHSAndEnergy(const fourdst::composition::Composition &comp, double T9, double rho) const {
+std::expected<gridfire::engine::StepDerivatives<double>, gridfire::engine::EngineStatus> PyEngine::calculateRHSAndEnergy(
+    const fourdst::composition::CompositionAbstract &comp,
+    double T9,
+    double rho
+) const {
     PYBIND11_OVERRIDE_PURE(
-        PYBIND11_TYPE(std::expected<gridfire::StepDerivatives<double>, gridfire::expectations::StaleEngineError>),
-        gridfire::Engine,
+        PYBIND11_TYPE(std::expected<gridfire::engine::StepDerivatives<double>, gridfire::engine::EngineStatus>),
+        gridfire::engine::Engine,
         calculateRHSAndEnergy,
         comp, T9, rho
     );
@@ -65,19 +69,29 @@ const std::vector<fourdst::atomic::Species>& PyDynamicEngine::getNetworkSpecies(
 
     py::pybind11_fail("Tried to call pure virtual function \"DynamicEngine::getNetworkSpecies\"");
 }
-std::expected<gridfire::StepDerivatives<double>, gridfire::expectations::StaleEngineError> PyDynamicEngine::calculateRHSAndEnergy(const fourdst::composition::Composition &comp, double T9, double rho) const {
+
+
+std::expected<gridfire::engine::StepDerivatives<double>, gridfire::engine::EngineStatus> PyDynamicEngine::calculateRHSAndEnergy(
+    const fourdst::composition::CompositionAbstract &comp,
+    double T9,
+    double rho
+) const {
     PYBIND11_OVERRIDE_PURE(
-        PYBIND11_TYPE(std::expected<gridfire::StepDerivatives<double>, gridfire::expectations::StaleEngineError>),
-        gridfire::Engine,
+        PYBIND11_TYPE(std::expected<gridfire::engine::StepDerivatives<double>, gridfire::engine::EngineStatus>),
+        gridfire::engine::DynamicEngine,
         calculateRHSAndEnergy,
         comp, T9, rho
     );
 }
 
-void PyDynamicEngine::generateJacobianMatrix(const fourdst::composition::Composition& comp, double T9, double rho) const {
+gridfire::engine::NetworkJacobian PyDynamicEngine::generateJacobianMatrix(
+    const fourdst::composition::CompositionAbstract& comp,
+    double T9,
+    double rho
+) const {
     PYBIND11_OVERRIDE_PURE(
-        void,
-        gridfire::DynamicEngine,
+        gridfire::engine::NetworkJacobian,
+        gridfire::engine::DynamicEngine,
         generateJacobianMatrix,
         comp,
         T9,
@@ -85,15 +99,15 @@ void PyDynamicEngine::generateJacobianMatrix(const fourdst::composition::Composi
     );
 }
 
-void PyDynamicEngine::generateJacobianMatrix(
-    const fourdst::composition::Composition &comp,
+gridfire::engine::NetworkJacobian PyDynamicEngine::generateJacobianMatrix(
+    const fourdst::composition::CompositionAbstract &comp,
     const double T9,
     const double rho,
     const std::vector<fourdst::atomic::Species> &activeSpecies
 ) const {
     PYBIND11_OVERRIDE_PURE(
-        void,
-        gridfire::DynamicEngine,
+        gridfire::engine::NetworkJacobian,
+        gridfire::engine::DynamicEngine,
         generateJacobianMatrix,
         comp,
         T9,
@@ -102,10 +116,15 @@ void PyDynamicEngine::generateJacobianMatrix(
     );
 }
 
-void PyDynamicEngine::generateJacobianMatrix(const fourdst::composition::Composition &comp, double T9, double rho, const gridfire::SparsityPattern &sparsityPattern) const {
+gridfire::engine::NetworkJacobian PyDynamicEngine::generateJacobianMatrix(
+    const fourdst::composition::CompositionAbstract &comp,
+    double T9,
+    double rho,
+    const gridfire::engine::SparsityPattern &sparsityPattern
+) const {
     PYBIND11_OVERRIDE_PURE(
-        void,
-        gridfire::DynamicEngine,
+        gridfire::engine::NetworkJacobian,
+        gridfire::engine::DynamicEngine,
         generateJacobianMatrix,
         comp,
         T9,
@@ -114,38 +133,36 @@ void PyDynamicEngine::generateJacobianMatrix(const fourdst::composition::Composi
     );
 }
 
-double PyDynamicEngine::getJacobianMatrixEntry(const fourdst::atomic::Species& rowSpecies, const fourdst::atomic::Species& colSpecies) const {
-    PYBIND11_OVERRIDE_PURE(
-        double,
-        gridfire::DynamicEngine,
-        getJacobianMatrixEntry,
-        rowSpecies,
-        colSpecies
-    );
-}
-
 void PyDynamicEngine::generateStoichiometryMatrix() {
     PYBIND11_OVERRIDE_PURE(
         void,
-        gridfire::DynamicEngine,
+        gridfire::engine::DynamicEngine,
         generateStoichiometryMatrix
     );
 }
 
-int PyDynamicEngine::getStoichiometryMatrixEntry(const fourdst::atomic::Species& species, const gridfire::reaction::Reaction& reaction) const {
+int PyDynamicEngine::getStoichiometryMatrixEntry(
+    const fourdst::atomic::Species& species,
+    const gridfire::reaction::Reaction& reaction
+) const {
     PYBIND11_OVERRIDE_PURE(
         int,
-        gridfire::DynamicEngine,
+        gridfire::engine::DynamicEngine,
         getStoichiometryMatrixEntry,
         species,
         reaction
     );
 }
 
-double PyDynamicEngine::calculateMolarReactionFlow(const gridfire::reaction::Reaction &reaction, const fourdst::composition::Composition &comp, double T9, double rho) const {
+double PyDynamicEngine::calculateMolarReactionFlow(
+    const gridfire::reaction::Reaction &reaction,
+    const fourdst::composition::CompositionAbstract &comp,
+    double T9,
+    double rho
+) const {
     PYBIND11_OVERRIDE_PURE(
         double,
-        gridfire::DynamicEngine,
+        gridfire::engine::DynamicEngine,
         calculateMolarReactionFlow,
         reaction,
         comp,
@@ -157,7 +174,7 @@ double PyDynamicEngine::calculateMolarReactionFlow(const gridfire::reaction::Rea
 const gridfire::reaction::ReactionSet& PyDynamicEngine::getNetworkReactions() const {
     PYBIND11_OVERRIDE_PURE(
         const gridfire::reaction::ReactionSet&,
-        gridfire::DynamicEngine,
+        gridfire::engine::DynamicEngine,
         getNetworkReactions
     );
 }
@@ -165,16 +182,20 @@ const gridfire::reaction::ReactionSet& PyDynamicEngine::getNetworkReactions() co
 void PyDynamicEngine::setNetworkReactions(const gridfire::reaction::ReactionSet& reactions) {
     PYBIND11_OVERRIDE_PURE(
         void,
-        gridfire::DynamicEngine,
+        gridfire::engine::DynamicEngine,
         setNetworkReactions,
         reactions
     );
 }
 
-std::expected<std::unordered_map<fourdst::atomic::Species, double>, gridfire::expectations::StaleEngineError> PyDynamicEngine::getSpeciesTimescales(const fourdst::composition::Composition &comp, double T9, double rho) const {
+std::expected<std::unordered_map<fourdst::atomic::Species, double>, gridfire::engine::EngineStatus> PyDynamicEngine::getSpeciesTimescales(
+    const fourdst::composition::CompositionAbstract &comp,
+    double T9,
+    double rho
+) const {
     PYBIND11_OVERRIDE_PURE(
-        PYBIND11_TYPE(std::expected<std::unordered_map<fourdst::atomic::Species, double>, gridfire::expectations::StaleEngineError>),
-        gridfire::DynamicEngine,
+        PYBIND11_TYPE(std::expected<std::unordered_map<fourdst::atomic::Species, double>, gridfire::engine::EngineStatus>),
+        gridfire::engine::DynamicEngine,
         getSpeciesTimescales,
         comp,
         T9,
@@ -182,10 +203,14 @@ std::expected<std::unordered_map<fourdst::atomic::Species, double>, gridfire::ex
     );
 }
 
-std::expected<std::unordered_map<fourdst::atomic::Species, double>, gridfire::expectations::StaleEngineError> PyDynamicEngine::getSpeciesDestructionTimescales(const fourdst::composition::Composition &comp, double T9, double rho) const {
+std::expected<std::unordered_map<fourdst::atomic::Species, double>, gridfire::engine::EngineStatus> PyDynamicEngine::getSpeciesDestructionTimescales(
+    const fourdst::composition::CompositionAbstract &comp,
+    double T9,
+    double rho
+) const {
     PYBIND11_OVERRIDE_PURE(
-        PYBIND11_TYPE(std::expected<std::unordered_map<fourdst::atomic::Species, double>, gridfire::expectations::StaleEngineError>),
-        gridfire::DynamicEngine,
+        PYBIND11_TYPE(std::expected<std::unordered_map<fourdst::atomic::Species, double>, gridfire::engine::EngineStatus>),
+        gridfire::engine::DynamicEngine,
         getSpeciesDestructionTimescales,
         comp, T9, rho
     );
@@ -194,7 +219,7 @@ std::expected<std::unordered_map<fourdst::atomic::Species, double>, gridfire::ex
 fourdst::composition::Composition PyDynamicEngine::update(const gridfire::NetIn &netIn) {
     PYBIND11_OVERRIDE_PURE(
         fourdst::composition::Composition,
-        gridfire::DynamicEngine,
+        gridfire::engine::DynamicEngine,
         update,
         netIn
     );
@@ -203,7 +228,7 @@ fourdst::composition::Composition PyDynamicEngine::update(const gridfire::NetIn
 bool PyDynamicEngine::isStale(const gridfire::NetIn &netIn) {
     PYBIND11_OVERRIDE_PURE(
         bool,
-        gridfire::DynamicEngine,
+        gridfire::engine::DynamicEngine,
         isStale,
         netIn
     );
@@ -212,7 +237,7 @@ bool PyDynamicEngine::isStale(const gridfire::NetIn &netIn) {
 void PyDynamicEngine::setScreeningModel(gridfire::screening::ScreeningType model) {
     PYBIND11_OVERRIDE_PURE(
         void,
-        gridfire::DynamicEngine,
+        gridfire::engine::DynamicEngine,
         setScreeningModel,
         model
     );
@@ -221,7 +246,7 @@ void PyDynamicEngine::setScreeningModel(gridfire::screening::ScreeningType model
 gridfire::screening::ScreeningType PyDynamicEngine::getScreeningModel() const {
     PYBIND11_OVERRIDE_PURE(
         gridfire::screening::ScreeningType,
-        gridfire::DynamicEngine,
+        gridfire::engine::DynamicEngine,
         getScreeningModel
     );
 }
@@ -229,7 +254,7 @@ gridfire::screening::ScreeningType PyDynamicEngine::getScreeningModel() const {
 size_t PyDynamicEngine::getSpeciesIndex(const fourdst::atomic::Species &species) const {
     PYBIND11_OVERRIDE_PURE(
         int,
-        gridfire::DynamicEngine,
+        gridfire::engine::DynamicEngine,
         getSpeciesIndex,
         species
     );
@@ -238,28 +263,28 @@ size_t PyDynamicEngine::getSpeciesIndex(const fourdst::atomic::Species &species)
 std::vector<double> PyDynamicEngine::mapNetInToMolarAbundanceVector(const gridfire::NetIn &netIn) const {
     PYBIND11_OVERRIDE_PURE(
         std::vector<double>,
-        gridfire::DynamicEngine,
+        gridfire::engine::DynamicEngine,
         mapNetInToMolarAbundanceVector,
         netIn
     );
 }
 
-gridfire::PrimingReport PyDynamicEngine::primeEngine(const gridfire::NetIn &netIn) {
+gridfire::engine::PrimingReport PyDynamicEngine::primeEngine(const gridfire::NetIn &netIn) {
     PYBIND11_OVERRIDE_PURE(
-        gridfire::PrimingReport,
-        gridfire::DynamicEngine,
+        gridfire::engine::PrimingReport,
+        gridfire::engine::DynamicEngine,
         primeEngine,
         netIn
     );
 }
 
-gridfire::EnergyDerivatives PyDynamicEngine::calculateEpsDerivatives(
-    const fourdst::composition::Composition &comp,
+gridfire::engine::EnergyDerivatives PyDynamicEngine::calculateEpsDerivatives(
+    const fourdst::composition::CompositionAbstract &comp,
     const double T9,
     const double rho) const {
     PYBIND11_OVERRIDE_PURE(
-        gridfire::EnergyDerivatives,
-        gridfire::DynamicEngine,
+        gridfire::engine::EnergyDerivatives,
+        gridfire::engine::DynamicEngine,
         calculateEpsDerivatives,
         comp,
         T9,
@@ -268,28 +293,41 @@ gridfire::EnergyDerivatives PyDynamicEngine::calculateEpsDerivatives(
 }
 
 fourdst::composition::Composition PyDynamicEngine::collectComposition(
-    fourdst::composition::Composition &comp
+    const fourdst::composition::CompositionAbstract &comp,
+    const double T9,
+    const double rho
 ) const {
     PYBIND11_OVERRIDE_PURE(
         fourdst::composition::Composition,
-        gridfire::DynamicEngine,
+        gridfire::engine::DynamicEngine,
         collectComposition,
-        comp
+        comp,
+        T9,
+        rho
     );
 }
 
-const gridfire::Engine& PyEngineView::getBaseEngine() const {
+gridfire::engine::SpeciesStatus PyDynamicEngine::getSpeciesStatus(const fourdst::atomic::Species &species) const {
     PYBIND11_OVERRIDE_PURE(
-        const gridfire::Engine&,
-        gridfire::EngineView<gridfire::Engine>,
+        gridfire::engine::SpeciesStatus,
+        gridfire::engine::DynamicEngine,
+        getSpeciesStatus,
+        species
+    );
+}
+
+const gridfire::engine::Engine& PyEngineView::getBaseEngine() const {
+    PYBIND11_OVERRIDE_PURE(
+        const gridfire::engine::Engine&,
+        gridfire::engine::EngineView<gridfire::engine::Engine>,
         getBaseEngine
     );
 }
 
-const gridfire::DynamicEngine& PyDynamicEngineView::getBaseEngine() const {
+const gridfire::engine::DynamicEngine& PyDynamicEngineView::getBaseEngine() const {
     PYBIND11_OVERRIDE_PURE(
-        const gridfire::DynamicEngine&,
-        gridfire::EngineView<gridfire::DynamicEngine>,
+        const gridfire::engine::DynamicEngine&,
+        gridfire::engine::EngineView<gridfire::engine::DynamicEngine>,
         getBaseEngine
     );
 }
diff --git a/src/python/engine/trampoline/py_engine.h b/src/python/engine/trampoline/py_engine.h
index e8f8d383..3eaaf3ff 100644
--- a/src/python/engine/trampoline/py_engine.h
+++ b/src/python/engine/trampoline/py_engine.h
@@ -1,7 +1,6 @@
 #pragma once
 
 #include "gridfire/engine/engine.h"
-#include "gridfire/expectations/expected_engine.h"
 
 #include "fourdst/atomic/atomicSpecies.h"
 
@@ -9,11 +8,11 @@
 #include <expected>
 
 
-class PyEngine final : public gridfire::Engine {
+class PyEngine final : public gridfire::engine::Engine {
 public:
     const std::vector<fourdst::atomic::Species>& getNetworkSpecies() const override;
 
-    std::expected<gridfire::StepDerivatives<double>,gridfire::expectations::StaleEngineError> calculateRHSAndEnergy(
+    std::expected<gridfire::engine::StepDerivatives<double>, gridfire::engine::EngineStatus> calculateRHSAndEnergy(
         const fourdst::composition::CompositionAbstract &comp,
         double T9,
         double rho
@@ -22,39 +21,34 @@ private:
     mutable std::vector<fourdst::atomic::Species> m_species_cache;
 };
 
-class PyDynamicEngine final : public gridfire::DynamicEngine {
+class PyDynamicEngine final : public gridfire::engine::DynamicEngine {
 public:
     const std::vector<fourdst::atomic::Species>& getNetworkSpecies() const override;
 
-    std::expected<gridfire::StepDerivatives<double>, gridfire::expectations::StaleEngineError> calculateRHSAndEnergy(
+    std::expected<gridfire::engine::StepDerivatives<double>, gridfire::engine::EngineStatus> calculateRHSAndEnergy(
         const fourdst::composition::CompositionAbstract &comp,
         double T9,
         double rho
     ) const override;
 
-    void generateJacobianMatrix(
-        const fourdst::composition::Composition& comp,
+    gridfire::engine::NetworkJacobian generateJacobianMatrix(
+        const fourdst::composition::CompositionAbstract& comp,
         double T9,
         double rho
     ) const override;
 
-    void generateJacobianMatrix(
+    gridfire::engine::NetworkJacobian generateJacobianMatrix(
         const fourdst::composition::CompositionAbstract &comp,
         double T9,
         double rho,
         const std::vector<fourdst::atomic::Species> &activeSpecies
     ) const override;
 
-    void generateJacobianMatrix(
-        const fourdst::composition::Composition& comp,
+    gridfire::engine::NetworkJacobian generateJacobianMatrix(
+        const fourdst::composition::CompositionAbstract& comp,
         double T9,
         double rho,
-        const gridfire::SparsityPattern &sparsityPattern
-    ) const override;
-
-    double getJacobianMatrixEntry(
-        const fourdst::atomic::Species& rowSpecies,
-        const fourdst::atomic::Species& colSpecies
+        const gridfire::engine::SparsityPattern &sparsityPattern
     ) const override;
 
     void generateStoichiometryMatrix() override;
@@ -77,14 +71,14 @@ public:
         const gridfire::reaction::ReactionSet& reactions
     ) override;
 
-    std::expected<std::unordered_map<fourdst::atomic::Species, double>, gridfire::expectations::StaleEngineError> getSpeciesTimescales(
+    std::expected<std::unordered_map<fourdst::atomic::Species, double>, gridfire::engine::EngineStatus> getSpeciesTimescales(
         const fourdst::composition::CompositionAbstract &comp,
         double T9,
         double rho
     ) const override;
 
-    std::expected<std::unordered_map<fourdst::atomic::Species, double>, gridfire::expectations::StaleEngineError> getSpeciesDestructionTimescales(
-        const fourdst::composition::Composition &comp,
+    std::expected<std::unordered_map<fourdst::atomic::Species, double>, gridfire::engine::EngineStatus> getSpeciesDestructionTimescales(
+        const fourdst::composition::CompositionAbstract &comp,
         double T9,
         double rho
     ) const override;
@@ -111,38 +105,44 @@ public:
         const gridfire::NetIn &netIn
     ) const override;
 
-    gridfire::PrimingReport primeEngine(
+    gridfire::engine::PrimingReport primeEngine(
         const gridfire::NetIn &netIn
     ) override;
 
-    gridfire::BuildDepthType getDepth() const override {
+    gridfire::engine::BuildDepthType getDepth() const override {
         throw std::logic_error("Network depth not supported by this engine.");
     }
     void rebuild(
         const fourdst::composition::CompositionAbstract &comp,
-        gridfire::BuildDepthType depth
+        gridfire::engine::BuildDepthType depth
     ) override {
         throw std::logic_error("Setting network depth not supported by this engine.");
     }
 
-    [[nodiscard]] gridfire::EnergyDerivatives calculateEpsDerivatives(
-        const fourdst::composition::Composition &comp,
+    [[nodiscard]] gridfire::engine::EnergyDerivatives calculateEpsDerivatives(
+        const fourdst::composition::CompositionAbstract &comp,
         double T9,
         double rho
     ) const override;
 
     fourdst::composition::Composition collectComposition(
-        fourdst::composition::Composition &comp
+        const fourdst::composition::CompositionAbstract &comp,
+        double T9,
+        double rho
+    ) const override;
+
+    gridfire::engine::SpeciesStatus getSpeciesStatus(
+        const fourdst::atomic::Species &species
     ) const override;
 
 private:
     mutable std::vector<fourdst::atomic::Species> m_species_cache;
 };
 
-class PyEngineView final : public gridfire::EngineView<gridfire::Engine> {
-    [[nodiscard]] const gridfire::Engine& getBaseEngine() const override;
+class PyEngineView final : public gridfire::engine::EngineView<gridfire::engine::Engine> {
+    [[nodiscard]] const gridfire::engine::Engine& getBaseEngine() const override;
 };
 
-class PyDynamicEngineView final : public gridfire::EngineView<gridfire::DynamicEngine> {
-    [[nodiscard]] const gridfire::DynamicEngine& getBaseEngine() const override;
+class PyDynamicEngineView final : public gridfire::engine::EngineView<gridfire::engine::DynamicEngine> {
+    [[nodiscard]] const gridfire::engine::DynamicEngine& getBaseEngine() const override;
 };
\ No newline at end of file
diff --git a/src/python/exceptions/bindings.cpp b/src/python/exceptions/bindings.cpp
index 3223fb0c..14110dce 100644
--- a/src/python/exceptions/bindings.cpp
+++ b/src/python/exceptions/bindings.cpp
@@ -7,48 +7,41 @@ namespace py = pybind11;
 #include "gridfire/exceptions/exceptions.h"
 
 void register_exception_bindings(const py::module &m) {
-    py::register_exception<gridfire::exceptions::EngineError>(m, "GridFireEngineError");
+    py::register_exception<gridfire::exceptions::GridFireError>(m, "GridFireError");
 
-    // TODO: Make it so that we can grab the stale state in python
-    // m.attr("StaleEngineTrigger") = py::register_exception<gridfire::exceptions::StaleEngineTrigger>(m, "StaleEngineTrigger", m.attr("GridFireEngineError"));
-    m.attr("StaleEngineError") = py::register_exception<gridfire::exceptions::StaleEngineError>(m, "StaleEngineError", m.attr("GridFireEngineError"));
-    m.attr("FailedToPartitionEngineError") = py::register_exception<gridfire::exceptions::FailedToPartitionEngineError>(m, "FailedToPartitionEngineError", m.attr("GridFireEngineError"));
-    m.attr("NetworkResizedError") = py::register_exception<gridfire::exceptions::NetworkResizedError>(m, "NetworkResizedError", m.attr("GridFireEngineError"));
-    m.attr("UnableToSetNetworkReactionsError") = py::register_exception<gridfire::exceptions::UnableToSetNetworkReactionsError>(m, "UnableToSetNetworkReactionsError", m.attr("GridFireEngineError"));
+    py::register_exception<gridfire::exceptions::DebugException>(m, "DebugException", m.attr("GridFireError"));
 
-    py::class_<gridfire::exceptions::StaleEngineTrigger::state>(m, "StaleEngineState")
-        .def(py::init<>())
-        .def_readwrite("T9", &gridfire::exceptions::StaleEngineTrigger::state::m_T9)
-        .def_readwrite("rho", &gridfire::exceptions::StaleEngineTrigger::state::m_rho)
-        .def_readwrite("Y", &gridfire::exceptions::StaleEngineTrigger::state::m_Y)
-        .def_readwrite("t", &gridfire::exceptions::StaleEngineTrigger::state::m_t)
-        .def_readwrite("total_steps", &gridfire::exceptions::StaleEngineTrigger::state::m_total_steps)
-        .def_readwrite("eps_nuc", &gridfire::exceptions::StaleEngineTrigger::state::m_eps_nuc);
+    py::register_exception<gridfire::exceptions::EngineError>(m, "EngineError", m.attr("GridFireError"));
 
-    py::class_<gridfire::exceptions::StaleEngineTrigger>(m, "StaleEngineTrigger")
-        .def(py::init<const gridfire::exceptions::StaleEngineTrigger::state &>())
-        .def("getState", &gridfire::exceptions::StaleEngineTrigger::getState)
-        .def("numSpecies", &gridfire::exceptions::StaleEngineTrigger::numSpecies)
-        .def("totalSteps", &gridfire::exceptions::StaleEngineTrigger::totalSteps)
-        .def("energy", &gridfire::exceptions::StaleEngineTrigger::energy)
-        .def("getMolarAbundance", &gridfire::exceptions::StaleEngineTrigger::getMolarAbundance)
-        .def("temperature", &gridfire::exceptions::StaleEngineTrigger::temperature)
-        .def("density", &gridfire::exceptions::StaleEngineTrigger::density)
-        .def("__repr__", [&](const gridfire::exceptions::StaleEngineTrigger& self) {
-           return self.what();
-        });
+    py::register_exception<gridfire::exceptions::FailedToPartitionEngineError>(m, "FailedToPartitionEngineError", m.attr("EngineError"));
+    py::register_exception<gridfire::exceptions::NetworkResizedError>(m, "NetworkResizedError", m.attr("EngineError"));
+    py::register_exception<gridfire::exceptions::UnableToSetNetworkReactionsError>(m, "UnableToSetNetworkReactionsError", m.attr("EngineError"));
+    py::register_exception<gridfire::exceptions::BadCollectionError>(m, "BadCollectionError", m.attr("EngineError"));
+    py::register_exception<gridfire::exceptions::InvalidQSESolutionError>(m, "InvalidQSESolutionError", m.attr("EngineError"));
+    py::register_exception<gridfire::exceptions::BadRHSEngineError>(m, "BadRHSEngineError", m.attr("EngineError"));
 
-    py::register_exception<gridfire::exceptions::FailedToPartitionEngineError>(m, "FailedToPartitionEngineError", m.attr("GridFireEngineError"));
-    py::register_exception<gridfire::exceptions::NetworkResizedError>(m, "NetworkResizedError", m.attr("GridFireEngineError"));
-    py::register_exception<gridfire::exceptions::UnableToSetNetworkReactionsError>(m, "UnableToSetNetworkReactionsError", m.attr("GridFireEngineError"));
-    py::register_exception<gridfire::exceptions::BadCollectionError>(m, "BadCollectionError", m.attr("GridFireEngineError"));
-
-    py::register_exception<gridfire::exceptions::JacobianError>(m, "JacobianError", m.attr("GridFireEngineError"));
+    py::register_exception<gridfire::exceptions::JacobianError>(m, "JacobianError", m.attr("EngineError"));
 
     py::register_exception<gridfire::exceptions::StaleJacobianError>(m, "StaleJacobianError", m.attr("JacobianError"));
     py::register_exception<gridfire::exceptions::UninitializedJacobianError>(m, "UninitializedJacobianError", m.attr("JacobianError"));
     py::register_exception<gridfire::exceptions::UnknownJacobianError>(m, "UnknownJacobianError", m.attr("JacobianError"));
 
-    py::register_exception<gridfire::exceptions::UtilityError>(m, "UtilityError");
+    py::register_exception<gridfire::exceptions::UtilityError>(m, "UtilityError", m.attr("GridFireError"));
     py::register_exception<gridfire::exceptions::HashingError>(m, "HashingError", m.attr("UtilityError"));
+
+    py::register_exception<gridfire::exceptions::PolicyError>(m, "PolicyError", m.attr("GridFireError"));
+    py::register_exception<gridfire::exceptions::MissingBaseReactionError>(m, "MissingBaseReactionError", m.attr("PolicyError"));
+    py::register_exception<gridfire::exceptions::MissingSeedSpeciesError>(m, "MissingSeedSpeciesError", m.attr("PolicyError"));
+    py::register_exception<gridfire::exceptions::MissingKeyReactionError>(m, "MissingKeyReactionError", m.attr("PolicyError"));
+
+    py::register_exception<gridfire::exceptions::ReactionError>(m, "ReactionError", m.attr("GridFireError"));
+    py::register_exception<gridfire::exceptions::ReactionParsingError>(m, "ReactionParsingError", m.attr("ReactionError"));
+
+    py::register_exception<gridfire::exceptions::SolverError>(m, "SolverError", m.attr("GridFireError"));
+    py::register_exception<gridfire::exceptions::SingularJacobianError>(m, "SingularJacobianError", m.attr("SolverError"));
+    py::register_exception<gridfire::exceptions::IllConditionedJacobianError>(m, "IllConditionedJacobianError", m.attr("SolverError"));
+    py::register_exception<gridfire::exceptions::SUNDIALSError>(m, "SUNDIALSError", m.attr("SolverError"));
+    py::register_exception<gridfire::exceptions::CVODESolverFailureError>(m, "CVODESolverFailureError", m.attr("SUNDIALSError"));
+    py::register_exception<gridfire::exceptions::KINSolSolverFailureError>(m, "KINSolSolverFailureError", m.attr("SUNDIALSError"));
+
 }
diff --git a/src/python/expectations/bindings.cpp b/src/python/expectations/bindings.cpp
deleted file mode 100644
index 7ff4c8bc..00000000
--- a/src/python/expectations/bindings.cpp
+++ /dev/null
@@ -1,43 +0,0 @@
-#include <pybind11/pybind11.h>
-#include <pybind11/stl.h> // Needed for vectors, maps, sets, strings
-#include <pybind11/stl_bind.h> // Needed for binding std::vector, std::map etc. if needed directly
-
-#include "bindings.h"
-
-namespace py = pybind11;
-
-#include "gridfire/expectations/expectations.h"
-
-void register_expectation_bindings(const py::module &m) {
-    py::enum_<gridfire::expectations::EngineErrorTypes>(m, "EngineErrorTypes")
-        .value("FAILURE", gridfire::expectations::EngineErrorTypes::FAILURE)
-        .value("INDEX", gridfire::expectations::EngineErrorTypes::INDEX)
-        .value("STALE", gridfire::expectations::EngineErrorTypes::STALE)
-        .export_values();
-
-    py::enum_<gridfire::expectations::StaleEngineErrorTypes>(m, "StaleEngineErrorTypes")
-        .value("SYSTEM_RESIZED", gridfire::expectations::StaleEngineErrorTypes::SYSTEM_RESIZED)
-        .export_values();
-
-    // Bind the base class
-    py::class_<gridfire::expectations::EngineError>(m, "EngineError")
-        .def_readonly("message", &gridfire::expectations::EngineError::m_message)
-        .def_readonly("type", &gridfire::expectations::EngineError::type)
-        .def("__str__", [](const gridfire::expectations::EngineError &e) {return e.m_message;});
-
-    // Bind the EngineIndexError, specifying EngineError as the base
-    py::class_<gridfire::expectations::EngineIndexError, gridfire::expectations::EngineError>(m, "EngineIndexError")
-        .def(py::init<int>(), py::arg("index"))
-        .def_readonly("index", &gridfire::expectations::EngineIndexError::m_index)
-        .def("__str__", [](const gridfire::expectations::EngineIndexError &e) {
-            return e.m_message + " at index " + std::to_string(e.m_index);
-        });
-
-    // Bind the StaleEngineError, specifying EngineError as the base
-    py::class_<gridfire::expectations::StaleEngineError, gridfire::expectations::EngineError>(m, "StaleEngineError")
-        .def(py::init<gridfire::expectations::StaleEngineErrorTypes>(), py::arg("stale_type"))
-        .def_readonly("stale_type", &gridfire::expectations::StaleEngineError::staleType)
-        .def("__str__", [](const gridfire::expectations::StaleEngineError &e) {
-            return static_cast<std::string>(e);
-        });
-}
diff --git a/src/python/expectations/bindings.h b/src/python/expectations/bindings.h
deleted file mode 100644
index f773583b..00000000
--- a/src/python/expectations/bindings.h
+++ /dev/null
@@ -1,5 +0,0 @@
-#pragma once
-
-#include <pybind11/pybind11.h>
-
-void register_expectation_bindings(const pybind11::module &m);
diff --git a/src/python/gridfire/__init__.py b/src/python/gridfire/__init__.py
new file mode 100644
index 00000000..d92f0f26
--- /dev/null
+++ b/src/python/gridfire/__init__.py
@@ -0,0 +1,20 @@
+from ._gridfire import *
+import sys
+
+from ._gridfire import type, utils, engine, solver, exceptions, partition, reaction, screening, io, policy
+
+sys.modules['gridfire.type'] = type
+sys.modules['gridfire.utils'] = utils
+sys.modules['gridfire.engine'] = engine
+sys.modules['gridfire.solver'] = solver
+sys.modules['gridfire.exceptions'] = exceptions
+sys.modules['gridfire.partition'] = partition
+sys.modules['gridfire.reaction'] = reaction
+sys.modules['gridfire.screening'] = screening
+sys.modules['gridfire.policy'] = policy
+sys.modules['gridfire.io'] = io
+
+__all__ = ['type', 'utils', 'engine', 'solver', 'exceptions', 'partition', 'reaction', 'screening', 'io', 'policy']
+
+__version__ = "v0.7.0_alpha_2025_10_25"
+
diff --git a/src/python/io/bindings.cpp b/src/python/io/bindings.cpp
index 21cb7bf0..00641787 100644
--- a/src/python/io/bindings.cpp
+++ b/src/python/io/bindings.cpp
@@ -19,9 +19,4 @@ auto register_io_bindings(const py::module &m) -> void {
         .def("parse", &gridfire::io::SimpleReactionListFileParser::parse,
             py::arg("filename"),
             "Parse a simple reaction list file and return a ParsedNetworkData object.");
-
-    // py::class_<gridfire::io::MESANetworkFileParser, gridfire::io::NetworkFileParser>(m, "MESANetworkFileParser")
-    //     .def("parse", &gridfire::io::MESANetworkFileParser::parse,
-    //         py::arg("filename"),
-    //         "Parse a MESA network file and return a ParsedNetworkData object.");
 }
\ No newline at end of file
diff --git a/src/python/meson.build b/src/python/meson.build
index ec84b6be..c38cfbed 100644
--- a/src/python/meson.build
+++ b/src/python/meson.build
@@ -1,10 +1,10 @@
 subdir('types')
 subdir('utils')
-subdir('expectations')
 subdir('exceptions')
 subdir('io')
 subdir('partition')
 subdir('reaction')
 subdir('screening')
 subdir('engine')
+subdir('policy')
 subdir('solver')
diff --git a/src/python/policy/bindings.cpp b/src/python/policy/bindings.cpp
new file mode 100644
index 00000000..175936a9
--- /dev/null
+++ b/src/python/policy/bindings.cpp
@@ -0,0 +1,233 @@
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <pybind11/stl_bind.h>
+
+#include <memory>
+#include "bindings.h"
+#include "trampoline/py_policy.h"
+
+#include "gridfire/policy/policy.h"
+
+PYBIND11_DECLARE_HOLDER_TYPE(T, std::unique_ptr<T>, true) // Declare unique_ptr as a holder type for pybind11
+
+namespace py = pybind11;
+
+namespace {
+
+    template <typename T>
+    concept IsReactionChainPolicy = std::is_base_of_v<gridfire::policy::ReactionChainPolicy, T>;
+
+    template <typename T>
+    concept IsNetworkPolicy = std::is_base_of_v<gridfire::policy::NetworkPolicy, T>;
+
+    template<IsReactionChainPolicy T, IsReactionChainPolicy BaseT>
+    void registerReactionChainPolicyDefs(py::class_<T, BaseT>& pyClass) {
+        pyClass.def(
+            "get_reactions",
+            &T::get_reactions,
+            "Get the ReactionSet representing this reaction chain."
+        )
+        .def(
+            "contains",
+            py::overload_cast<const std::string&>(&T::contains, py::const_),
+            py::arg("id"),
+            "Check if the reaction chain contains a reaction with the given ID."
+        )
+        .def(
+            "contains",
+            py::overload_cast<const gridfire::reaction::Reaction&>(&T::contains, py::const_),
+            py::arg("reaction"),
+            "Check if the reaction chain contains the given reaction."
+        )
+        .def(
+            "name",
+            &T::name,
+            "Get the name of the reaction chain policy."
+        )
+        .def(
+            "hash",
+            &T::hash,
+            py::arg("seed"),
+            "Compute a hash value for the reaction chain policy."
+        )
+        .def(
+            "__eq__",
+            &T::operator==,
+            py::arg("other"),
+            "Check equality with another ReactionChainPolicy."
+        )
+        .def(
+            "__ne__",
+            &T::operator!=,
+            py::arg("other"),
+            "Check inequality with another ReactionChainPolicy."
+        )
+        .def("__hash__", [](const T &self) {
+            return self.hash(0);
+            }
+        )
+        .def("__repr__", [](const T &self) {
+            std::stringstream ss;
+            ss << self;
+            return ss.str();
+        });
+    }
+
+    template<IsNetworkPolicy T, IsNetworkPolicy BaseT>
+    void registerNetworkPolicyDefs(py::class_<T, BaseT> pyClass) {
+        pyClass.def(
+            "name",
+            &T::name,
+            "Get the name of the network policy."
+        )
+        .def(
+            "get_seed_species",
+            &T::get_seed_species,
+            "Get the set of seed species required by the network policy."
+        )
+        .def(
+            "get_seed_reactions",
+            &T::get_seed_reactions,
+            "Get the set of seed reactions required by the network policy."
+        )
+        .def(
+            "get_status",
+            &T::get_status,
+            "Get the current status of the network policy."
+        )
+        .def(
+            "get_engine_types_stack",
+            &T::get_engine_types_stack,
+            "Get the types of engines in the stack constructed by the network policy."
+        )
+        .def(
+            "construct",
+            &T::construct,
+            py::return_value_policy::reference,
+            "Construct the network according to the policy."
+
+        );
+    }
+}
+
+
+void register_policy_bindings(pybind11::module &m) {
+    register_reaction_chain_policy_bindings(m);
+    register_network_policy_bindings(m);
+}
+
+void register_reaction_chain_policy_bindings(pybind11::module &m) {
+    using namespace gridfire::policy;
+
+    py::class_<ReactionChainPolicy, PyReactionChainPolicy> py_reactionChainPolicy(m, "ReactionChainPolicy");
+    py::class_<MultiReactionChainPolicy, ReactionChainPolicy> py_multiChainPolicy(m, "MultiReactionChainPolicy");
+    py::class_<TemperatureDependentChainPolicy, ReactionChainPolicy> py_tempDepChainPolicy(m, "TemperatureDependentChainPolicy");
+
+
+    py::class_<ProtonProtonIChainPolicy, TemperatureDependentChainPolicy> py_ppI(m, "ProtonProtonIChainPolicy");
+    py_ppI.def(py::init<>());
+    py_ppI.def("name", &ProtonProtonIChainPolicy::name, "Get the name of the reaction chain policy.");
+
+    py::class_<ProtonProtonIIChainPolicy, TemperatureDependentChainPolicy> py_ppII(m, "ProtonProtonIIChainPolicy");
+    py_ppII.def(py::init<>());
+    py_ppII.def("name", &ProtonProtonIIChainPolicy::name, "Get the name of the reaction chain policy.");
+
+    py::class_<ProtonProtonIIIChainPolicy, TemperatureDependentChainPolicy> py_ppIII(m, "ProtonProtonIIIChainPolicy");
+    py_ppIII.def(py::init<>());
+    py_ppIII.def("name", &ProtonProtonIIIChainPolicy::name, "Get the name of the reaction chain policy.");
+
+    py::class_<ProtonProtonChainPolicy, MultiReactionChainPolicy> py_ppChain(m, "ProtonProtonChainPolicy");
+    py_ppChain.def(py::init<>());
+    py_ppChain.def("name", &ProtonProtonChainPolicy::name, "Get the name of the reaction chain policy.");
+
+    registerReactionChainPolicyDefs(py_ppI);
+    registerReactionChainPolicyDefs(py_ppII);
+    registerReactionChainPolicyDefs(py_ppIII);
+    registerReactionChainPolicyDefs(py_ppChain);
+
+    py::class_<CNOIChainPolicy, TemperatureDependentChainPolicy> py_cnoI(m, "CNOIChainPolicy");
+    py_cnoI.def(py::init<>());
+    py_cnoI.def("name", &CNOIChainPolicy::name, "Get the name of the reaction chain policy.");
+
+    py::class_<CNOIIChainPolicy, TemperatureDependentChainPolicy> py_cnoII(m, "CNOIIChainPolicy");
+    py_cnoII.def(py::init<>());
+    py_cnoII.def("name", &CNOIIChainPolicy::name, "Get the name of the reaction chain policy.");
+
+    py::class_<CNOIIIChainPolicy, TemperatureDependentChainPolicy> py_cnoIII(m, "CNOIIIChainPolicy");
+    py_cnoIII.def(py::init<>());
+    py_cnoIII.def("name", &CNOIIIChainPolicy::name, "Get the name of the reaction chain policy.");
+
+    py::class_<CNOIVChainPolicy, TemperatureDependentChainPolicy> py_cnoIV(m, "CNOIVChainPolicy");
+    py_cnoIV.def(py::init<>());
+    py_cnoIV.def("name", &CNOIVChainPolicy::name, "Get the name of the reaction chain policy.");
+
+    py::class_<CNOChainPolicy, MultiReactionChainPolicy> py_cnoChain(m, "CNOChainPolicy");
+    py_cnoChain.def(py::init<>());
+    py_cnoChain.def("name", &CNOChainPolicy::name, "Get the name of the reaction chain policy.");
+
+    registerReactionChainPolicyDefs(py_cnoI);
+    registerReactionChainPolicyDefs(py_cnoII);
+    registerReactionChainPolicyDefs(py_cnoIII);
+    registerReactionChainPolicyDefs(py_cnoIV);
+    registerReactionChainPolicyDefs(py_cnoChain);
+
+    py::class_<HotCNOIChainPolicy, TemperatureDependentChainPolicy> py_hotCNOI(m, "HotCNOIChainPolicy");
+    py_hotCNOI.def(py::init<>());
+    py_hotCNOI.def("name", &HotCNOIChainPolicy::name, "Get the name of the reaction chain policy.");
+
+    py::class_<HotCNOIIChainPolicy, TemperatureDependentChainPolicy> py_hotCNOII(m, "HotCNOIIChainPolicy");
+    py_hotCNOII.def(py::init<>());
+    py_hotCNOII.def("name", &HotCNOIIChainPolicy::name, "Get the name of the reaction chain policy.");
+
+    py::class_<HotCNOIIIChainPolicy, TemperatureDependentChainPolicy> py_hotCNOIII(m, "HotCNOIIIChainPolicy");
+    py_hotCNOIII.def(py::init<>());
+    py_hotCNOIII.def("name", &HotCNOIIIChainPolicy::name, "Get the name of the reaction chain policy.");
+
+    py::class_<HotCNOChainPolicy, MultiReactionChainPolicy> py_hotCNOChain(m, "HotCNOChainPolicy");
+    py_hotCNOChain.def(py::init<>());
+    py_hotCNOChain.def("name", &HotCNOChainPolicy::name, "Get the name of the reaction chain policy.");
+
+    registerReactionChainPolicyDefs(py_hotCNOI);
+    registerReactionChainPolicyDefs(py_hotCNOII);
+    registerReactionChainPolicyDefs(py_hotCNOIII);
+    registerReactionChainPolicyDefs(py_hotCNOChain);
+
+    py::class_<TripleAlphaChainPolicy, TemperatureDependentChainPolicy> py_tripleAlpha(m, "TripleAlphaChainPolicy");
+    py_tripleAlpha.def(py::init<>());
+    py_tripleAlpha.def("name", &TripleAlphaChainPolicy::name, "Get the name of the reaction chain policy.");
+
+    registerReactionChainPolicyDefs(py_tripleAlpha);
+
+    py::class_<MainSequenceReactionChainPolicy, MultiReactionChainPolicy> py_mainSeq(m, "MainSequenceReactionChainPolicy");
+    py_mainSeq.def(py::init<>());
+    py_mainSeq.def("name", &MainSequenceReactionChainPolicy::name, "Get the name of the reaction chain policy.");
+
+    registerReactionChainPolicyDefs(py_mainSeq);
+
+}
+
+void register_network_policy_bindings(pybind11::module &m) {
+    py::enum_<gridfire::policy::NetworkPolicyStatus>(m, "NetworkPolicyStatus")
+        .value("UNINITIALIZED", gridfire::policy::NetworkPolicyStatus::UNINITIALIZED)
+        .value("INITIALIZED_UNVERIFIED", gridfire::policy::NetworkPolicyStatus::INITIALIZED_UNVERIFIED)
+        .value("MISSING_KEY_REACTION", gridfire::policy::NetworkPolicyStatus::MISSING_KEY_REACTION)
+        .value("MISSING_KEY_SPECIES", gridfire::policy::NetworkPolicyStatus::MISSING_KEY_SPECIES)
+        .value("INITIALIZED_VERIFIED", gridfire::policy::NetworkPolicyStatus::INITIALIZED_VERIFIED)
+        .export_values();
+
+    py::class_<gridfire::policy::NetworkPolicy, PyNetworkPolicy> py_networkPolicy(m, "NetworkPolicy");
+    py::class_<gridfire::policy::MainSequencePolicy, gridfire::policy::NetworkPolicy> py_mainSeqPolicy(m, "MainSequencePolicy");
+    py_mainSeqPolicy.def(
+        py::init<const fourdst::composition::Composition&>(),
+        py::arg("composition"),
+        "Construct MainSequencePolicy from an existing composition."
+    );
+    py_mainSeqPolicy.def(
+        py::init<std::vector<fourdst::atomic::Species>, const std::vector<double>&>(),
+        py::arg("seed_species"),
+        py::arg("mass_fractions"),
+        "Construct MainSequencePolicy from seed species and mass fractions."
+    );
+
+    registerNetworkPolicyDefs(py_mainSeqPolicy);
+}
diff --git a/src/python/policy/bindings.h b/src/python/policy/bindings.h
new file mode 100644
index 00000000..64903302
--- /dev/null
+++ b/src/python/policy/bindings.h
@@ -0,0 +1,7 @@
+#pragma once
+
+#include <pybind11/pybind11.h>
+
+void register_policy_bindings(pybind11::module& m);
+void register_reaction_chain_policy_bindings(pybind11::module& m);
+void register_network_policy_bindings(pybind11::module& m);
\ No newline at end of file
diff --git a/src/python/expectations/meson.build b/src/python/policy/meson.build
similarity index 88%
rename from src/python/expectations/meson.build
rename to src/python/policy/meson.build
index 569cdf79..9913850e 100644
--- a/src/python/expectations/meson.build
+++ b/src/python/policy/meson.build
@@ -1,4 +1,6 @@
 # Define the library
+subdir('trampoline')
+
 bindings_sources = files('bindings.cpp')
 bindings_headers = files('bindings.h')
 
@@ -8,7 +10,7 @@ dependencies = [
     pybind11_dep,
 ]
 
-shared_module('py_gf_expectations',
+shared_module('py_gf_policy',
               bindings_sources,
               cpp_args: ['-fvisibility=default'],
               install : true,
diff --git a/src/python/policy/trampoline/meson.build b/src/python/policy/trampoline/meson.build
new file mode 100644
index 00000000..e0cf7622
--- /dev/null
+++ b/src/python/policy/trampoline/meson.build
@@ -0,0 +1,21 @@
+gf_policy_trampoline_sources = files('py_policy.cpp')
+
+gf_policy_trapoline_dependencies = [
+    gridfire_dep,
+    pybind11_dep,
+    python3_dep,
+]
+
+gf_policy_trampoline_lib = static_library(
+    'policy_trampolines',
+    gf_policy_trampoline_sources,
+    include_directories: include_directories('.'),
+    dependencies: gf_policy_trapoline_dependencies,
+    install: false,
+)
+
+gr_policy_trampoline_dep = declare_dependency(
+    link_with: gf_policy_trampoline_lib,
+    include_directories: ('.'),
+    dependencies: gf_policy_trapoline_dependencies,
+)
diff --git a/src/python/policy/trampoline/py_policy.cpp b/src/python/policy/trampoline/py_policy.cpp
new file mode 100644
index 00000000..600f4186
--- /dev/null
+++ b/src/python/policy/trampoline/py_policy.cpp
@@ -0,0 +1,149 @@
+#include "py_policy.h"
+
+#include "pybind11/pybind11.h"
+#include "pybind11/stl.h"
+
+#include "fourdst/atomic/atomicSpecies.h"
+
+#include "gridfire/reaction/reaction.h"
+#include "gridfire/engine/engine.h"
+
+#include "gridfire/policy/policy.h"
+
+#include <string>
+#include <set>
+
+namespace py = pybind11;
+
+std::string PyNetworkPolicy::name() const {
+    PYBIND11_OVERRIDE_PURE(
+        std::string,
+        gridfire::policy::NetworkPolicy,
+        name
+    );
+}
+
+const std::set<fourdst::atomic::Species>& PyNetworkPolicy::get_seed_species() const {
+    PYBIND11_OVERRIDE_PURE(
+        const std::set<fourdst::atomic::Species>&,
+        gridfire::policy::NetworkPolicy,
+        get_seed_species
+    );
+}
+
+const gridfire::reaction::ReactionSet& PyNetworkPolicy::get_seed_reactions() const {
+    PYBIND11_OVERRIDE_PURE(
+        const gridfire::reaction::ReactionSet&,
+        gridfire::policy::NetworkPolicy,
+        get_seed_reactions
+    );
+}
+
+gridfire::engine::DynamicEngine& PyNetworkPolicy::construct() {
+    PYBIND11_OVERRIDE_PURE(
+        gridfire::engine::DynamicEngine&,
+        gridfire::policy::NetworkPolicy,
+        construct
+    );
+}
+
+gridfire::policy::NetworkPolicyStatus PyNetworkPolicy::get_status() const {
+    PYBIND11_OVERRIDE_PURE(
+        gridfire::policy::NetworkPolicyStatus,
+        gridfire::policy::NetworkPolicy,
+        getStatus
+    );
+}
+
+const std::vector<std::unique_ptr<gridfire::engine::DynamicEngine>> &PyNetworkPolicy::get_engine_stack() const {
+    PYBIND11_OVERRIDE_PURE(
+        const std::vector<std::unique_ptr<gridfire::engine::DynamicEngine>> &,
+        gridfire::policy::NetworkPolicy,
+        get_engine_stack
+    );
+}
+
+std::vector<gridfire::engine::EngineTypes> PyNetworkPolicy::get_engine_types_stack() const {
+    PYBIND11_OVERRIDE_PURE(
+        std::vector<gridfire::engine::EngineTypes>,
+        gridfire::policy::NetworkPolicy,
+        get_engine_types_stack
+    );
+}
+
+const std::unique_ptr<gridfire::partition::PartitionFunction>& PyNetworkPolicy::get_partition_function() const {
+    PYBIND11_OVERRIDE_PURE(
+        const std::unique_ptr<gridfire::partition::PartitionFunction>&,
+        gridfire::policy::NetworkPolicy,
+        get_partition_function
+    );
+}
+
+const gridfire::reaction::ReactionSet &PyReactionChainPolicy::get_reactions() const {
+    PYBIND11_OVERRIDE_PURE(
+        const gridfire::reaction::ReactionSet &,
+        gridfire::policy::ReactionChainPolicy,
+        get_reactions
+    );
+}
+
+bool PyReactionChainPolicy::contains(const std::string &id) const {
+    PYBIND11_OVERRIDE_PURE(
+        bool,
+        gridfire::policy::ReactionChainPolicy,
+        contains,
+        id
+    );
+}
+
+bool PyReactionChainPolicy::contains(const gridfire::reaction::Reaction &reaction) const {
+    PYBIND11_OVERRIDE_PURE(
+        bool,
+        gridfire::policy::ReactionChainPolicy,
+        contains,
+        reaction
+    );
+}
+
+std::unique_ptr<gridfire::policy::ReactionChainPolicy> PyReactionChainPolicy::clone() const {
+    PYBIND11_OVERRIDE_PURE(
+        std::unique_ptr<gridfire::policy::ReactionChainPolicy>,
+        gridfire::policy::ReactionChainPolicy,
+        clone
+    );
+}
+
+std::string PyReactionChainPolicy::name() const {
+    PYBIND11_OVERRIDE_PURE(
+        std::string,
+        gridfire::policy::ReactionChainPolicy,
+        name
+    );
+}
+
+uint64_t PyReactionChainPolicy::hash(uint64_t seed) const {
+    PYBIND11_OVERRIDE_PURE(
+        uint64_t,
+        gridfire::policy::ReactionChainPolicy,
+        hash,
+        seed
+    );
+}
+
+bool PyReactionChainPolicy::operator==(const ReactionChainPolicy &other) const {
+    PYBIND11_OVERRIDE_PURE(
+        bool,
+        gridfire::policy::ReactionChainPolicy,
+        operator==,
+        other
+    );
+}
+
+bool PyReactionChainPolicy::operator!=(const ReactionChainPolicy &other) const {
+    PYBIND11_OVERRIDE_PURE(
+        bool,
+        gridfire::policy::ReactionChainPolicy,
+        operator!=,
+        other
+    );
+}
diff --git a/src/python/policy/trampoline/py_policy.h b/src/python/policy/trampoline/py_policy.h
new file mode 100644
index 00000000..91c6692b
--- /dev/null
+++ b/src/python/policy/trampoline/py_policy.h
@@ -0,0 +1,44 @@
+#pragma once
+
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+#include "gridfire/policy/policy.h"
+
+class PyNetworkPolicy final : public gridfire::policy::NetworkPolicy {
+public:
+    [[nodiscard]] std::string name() const override;
+
+    [[nodiscard]] const std::set<fourdst::atomic::Species>& get_seed_species() const override;
+
+    [[nodiscard]] const gridfire::reaction::ReactionSet& get_seed_reactions() const override;
+
+    [[nodiscard]] gridfire::engine::DynamicEngine& construct() override;
+
+    [[nodiscard]] gridfire::policy::NetworkPolicyStatus get_status() const override;
+
+    [[nodiscard]] const std::vector<std::unique_ptr<gridfire::engine::DynamicEngine>> &get_engine_stack() const override;
+
+    [[nodiscard]] std::vector<gridfire::engine::EngineTypes> get_engine_types_stack() const override;
+
+    [[nodiscard]] const std::unique_ptr<gridfire::partition::PartitionFunction>& get_partition_function() const override;
+};
+
+class PyReactionChainPolicy final : public gridfire::policy::ReactionChainPolicy {
+public:
+    [[nodiscard]] const gridfire::reaction::ReactionSet & get_reactions() const override;
+
+    [[nodiscard]] bool contains(const std::string &id) const override;
+
+    [[nodiscard]] bool contains(const gridfire::reaction::Reaction &reaction) const override;
+
+    [[nodiscard]] std::unique_ptr<ReactionChainPolicy> clone() const override;
+
+    [[nodiscard]] std::string name() const override;
+
+    [[nodiscard]] uint64_t hash(uint64_t seed) const override;
+
+    [[nodiscard]] bool operator==(const ReactionChainPolicy &other) const override;
+
+    [[nodiscard]] bool operator!=(const ReactionChainPolicy &other) const override;
+};
\ No newline at end of file
diff --git a/src/python/solver/bindings.cpp b/src/python/solver/bindings.cpp
index b00307b8..dd4e85f3 100644
--- a/src/python/solver/bindings.cpp
+++ b/src/python/solver/bindings.cpp
@@ -1,12 +1,10 @@
 #include <pybind11/pybind11.h>
-#include <pybind11/functional.h> // needed for std::function
 #include <pybind11/stl.h> // Needed for vectors, maps, sets, strings
 #include <pybind11/stl_bind.h> // Needed for binding std::vector, std::map etc. if needed directly
 #include <pybind11/numpy.h>
+#include <pybind11/functional.h>
 #include <functional>
 
-#include <boost/numeric/ublas/vector.hpp>
-
 #include "bindings.h"
 
 #include "gridfire/solver/strategies/CVODE_solver_strategy.h"
@@ -37,7 +35,7 @@ void register_solver_bindings(const py::module &m) {
     py_cvode_timestep_context.def_readonly("currentNonlinearIterations", &gridfire::solver::CVODESolverStrategy::TimestepContext::currentNonlinearIterations);
     py_cvode_timestep_context.def_property_readonly(
         "engine",
-        [](const gridfire::solver::CVODESolverStrategy::TimestepContext& self) -> const gridfire::DynamicEngine& {
+        [](const gridfire::solver::CVODESolverStrategy::TimestepContext& self) -> const gridfire::engine::DynamicEngine& {
             return self.engine;
         }
     );
@@ -66,7 +64,7 @@ void register_solver_bindings(const py::module &m) {
     auto py_cvode_solver_strategy = py::class_<gridfire::solver::CVODESolverStrategy, gridfire::solver::DynamicNetworkSolverStrategy>(m, "CVODESolverStrategy");
 
     py_cvode_solver_strategy.def(
-        py::init<gridfire::DynamicEngine&>(),
+        py::init<gridfire::engine::DynamicEngine&>(),
         py::arg("engine"),
         "Initialize the CVODESolverStrategy object."
     );
@@ -75,7 +73,7 @@ void register_solver_bindings(const py::module &m) {
         "evaluate",
         py::overload_cast<const gridfire::NetIn&, bool>(&gridfire::solver::CVODESolverStrategy::evaluate),
         py::arg("netIn"),
-        py::arg("display_trigger"),
+        py::arg("display_trigger") = false,
         "evaluate the dynamic engine using the dynamic engine class"
     );
 
@@ -92,6 +90,32 @@ void register_solver_bindings(const py::module &m) {
         "Enable logging to standard output."
     );
 
+    py_cvode_solver_strategy.def(
+        "set_absTol",
+        &gridfire::solver::CVODESolverStrategy::set_absTol,
+        py::arg("absTol"),
+        "Set the absolute tolerance for the CVODE solver."
+    );
+
+    py_cvode_solver_strategy.def(
+        "set_relTol",
+        &gridfire::solver::CVODESolverStrategy::set_relTol,
+        py::arg("relTol"),
+        "Set the relative tolerance for the CVODE solver."
+    );
+
+    py_cvode_solver_strategy.def(
+        "get_absTol",
+        &gridfire::solver::CVODESolverStrategy::get_absTol,
+        "Get the absolute tolerance for the CVODE solver."
+    );
+
+    py_cvode_solver_strategy.def(
+        "get_relTol",
+        &gridfire::solver::CVODESolverStrategy::get_relTol,
+        "Get the relative tolerance for the CVODE solver."
+    );
+
     py_cvode_solver_strategy.def(
     "set_callback",
     [](
@@ -103,7 +127,5 @@ void register_solver_bindings(const py::module &m) {
     py::arg("cb"),
     "Set a callback function which will run at the end of every successful timestep"
     );
-
-
 }
 
diff --git a/src/python/solver/trampoline/py_solver.h b/src/python/solver/trampoline/py_solver.h
index 25b19609..ef0713de 100644
--- a/src/python/solver/trampoline/py_solver.h
+++ b/src/python/solver/trampoline/py_solver.h
@@ -8,7 +8,7 @@
 #include <any>
 
 class PyDynamicNetworkSolverStrategy final : public gridfire::solver::DynamicNetworkSolverStrategy {
-    explicit PyDynamicNetworkSolverStrategy(gridfire::DynamicEngine &engine) : gridfire::solver::DynamicNetworkSolverStrategy(engine) {}
+    explicit PyDynamicNetworkSolverStrategy(gridfire::engine::DynamicEngine &engine) : gridfire::solver::DynamicNetworkSolverStrategy(engine) {}
     gridfire::NetOut evaluate(const gridfire::NetIn &netIn) override;
     void set_callback(const std::any &callback) override;
     [[nodiscard]] std::vector<std::tuple<std::string, std::string>> describe_callback_context() const override;
diff --git a/src/python/types/bindings.cpp b/src/python/types/bindings.cpp
index 23596f33..5e611c69 100644
--- a/src/python/types/bindings.cpp
+++ b/src/python/types/bindings.cpp
@@ -7,7 +7,7 @@
 
 namespace py = pybind11;
 
-#include "gridfire/network.h"
+#include "gridfire/types/types.h"
 
 void register_type_bindings(const pybind11::module &m) {
     py::class_<gridfire::NetIn>(m, "NetIn")
diff --git a/stubs/gridfire/__init__.pyi b/stubs/gridfire/__init__.pyi
new file mode 100644
index 00000000..e15f3ac9
--- /dev/null
+++ b/stubs/gridfire/__init__.pyi
@@ -0,0 +1,15 @@
+from __future__ import annotations
+from gridfire._gridfire import engine
+from gridfire._gridfire import exceptions
+from gridfire._gridfire import io
+from gridfire._gridfire import partition
+from gridfire._gridfire import policy
+from gridfire._gridfire import reaction
+from gridfire._gridfire import screening
+from gridfire._gridfire import solver
+from gridfire._gridfire import type
+from gridfire._gridfire import utils
+import sys as sys
+from . import _gridfire
+__all__: list = ['type', 'utils', 'engine', 'solver', 'exceptions', 'partition', 'reaction', 'screening', 'io', 'core', 'cli']
+__version__: str = 'v0.7.0_alpha_2025_10_25'
diff --git a/stubs/gridfire/_gridfire/__init__.pyi b/stubs/gridfire/_gridfire/__init__.pyi
new file mode 100644
index 00000000..23b26003
--- /dev/null
+++ b/stubs/gridfire/_gridfire/__init__.pyi
@@ -0,0 +1,15 @@
+"""
+Python bindings for the fourdst utility modules which are a part of the 4D-STAR project.
+"""
+from __future__ import annotations
+from . import engine
+from . import exceptions
+from . import io
+from . import partition
+from . import policy
+from . import reaction
+from . import screening
+from . import solver
+from . import type
+from . import utils
+__all__: list[str] = ['engine', 'exceptions', 'io', 'partition', 'policy', 'reaction', 'screening', 'solver', 'type', 'utils']
diff --git a/stubs/gridfire/_gridfire/engine/__init__.pyi b/stubs/gridfire/_gridfire/engine/__init__.pyi
new file mode 100644
index 00000000..b652bc70
--- /dev/null
+++ b/stubs/gridfire/_gridfire/engine/__init__.pyi
@@ -0,0 +1,1170 @@
+"""
+Engine and Engine View bindings
+"""
+from __future__ import annotations
+import collections.abc
+import fourdst._phys.atomic
+import fourdst._phys.composition
+import gridfire._gridfire.io
+import gridfire._gridfire.partition
+import gridfire._gridfire.reaction
+import gridfire._gridfire.screening
+import gridfire._gridfire.type
+import numpy
+import numpy.typing
+import typing
+from . import diagnostics
+__all__: list[str] = ['ACTIVE', 'ADAPTIVE_ENGINE_VIEW', 'AdaptiveEngineView', 'BuildDepthType', 'DEFAULT', 'DEFINED_ENGINE_VIEW', 'DefinedEngineView', 'DynamicEngine', 'EQUILIBRIUM', 'Engine', 'EngineTypes', 'FILE_DEFINED_ENGINE_VIEW', 'FULL_SUCCESS', 'FifthOrder', 'FileDefinedEngineView', 'FourthOrder', 'Full', 'GRAPH_ENGINE', 'GraphEngine', 'INACTIVE_FLOW', 'MAX_ITERATIONS_REACHED', 'MULTISCALE_PARTITIONING_ENGINE_VIEW', 'MultiscalePartitioningEngineView', 'NONE', 'NOT_PRESENT', 'NO_SPECIES_TO_PRIME', 'NetworkBuildDepth', 'NetworkConstructionFlags', 'NetworkJacobian', 'NetworkPrimingEngineView', 'PRIMING_ENGINE_VIEW', 'PrimingReport', 'PrimingReportStatus', 'REACLIB', 'REACLIB_STRONG', 'REACLIB_WEAK', 'SecondOrder', 'Shallow', 'SparsityPattern', 'SpeciesStatus', 'StepDerivatives', 'ThirdOrder', 'WRL_BETA_MINUS', 'WRL_BETA_PLUS', 'WRL_ELECTRON_CAPTURE', 'WRL_POSITRON_CAPTURE', 'WRL_WEAK', 'build_nuclear_network', 'diagnostics', 'primeNetwork', 'regularize_jacobian']
+class AdaptiveEngineView(DynamicEngine):
+    def __init__(self, baseEngine: DynamicEngine) -> None:
+        """
+        Construct an adaptive engine view with a base engine.
+        """
+    def calculateEpsDerivatives(self, comp: ..., T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> ...:
+        """
+        Calculate deps/dT and deps/drho
+        """
+    def calculateMolarReactionFlow(self: DynamicEngine, reaction: ..., comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> float:
+        """
+        Calculate the molar reaction flow for a given reaction.
+        """
+    def calculateRHSAndEnergy(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> StepDerivatives:
+        """
+        Calculate the right-hand side (dY/dt) and energy generation rate.
+        """
+    def collectComposition(self, composition: ..., T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> fourdst._phys.composition.Composition:
+        """
+        Recursively collect composition from current engine and any sub engines if they exist.
+        """
+    @typing.overload
+    def generateJacobianMatrix(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> NetworkJacobian:
+        """
+        Generate the Jacobian matrix for the current state.
+        """
+    @typing.overload
+    def generateJacobianMatrix(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat, activeSpecies: collections.abc.Sequence[fourdst._phys.atomic.Species]) -> NetworkJacobian:
+        """
+        Generate the jacobian matrix only for the subset of the matrix representing the active species.
+        """
+    @typing.overload
+    def generateJacobianMatrix(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat, sparsityPattern: collections.abc.Sequence[tuple[typing.SupportsInt, typing.SupportsInt]]) -> NetworkJacobian:
+        """
+        Generate the jacobian matrix for the given sparsity pattern
+        """
+    def generateStoichiometryMatrix(self) -> None:
+        ...
+    def getBaseEngine(self) -> DynamicEngine:
+        """
+        Get the base engine associated with this adaptive engine view.
+        """
+    def getDepth(self) -> gridfire._gridfire.engine.NetworkBuildDepth | int:
+        """
+        Get the current build depth of the engine.
+        """
+    def getNetworkReactions(self) -> gridfire._gridfire.reaction.ReactionSet:
+        """
+        Get the set of logical reactions in the network.
+        """
+    def getNetworkSpecies(self) -> list[fourdst._phys.atomic.Species]:
+        """
+        Get the list of species in the network.
+        """
+    def getScreeningModel(self) -> gridfire._gridfire.screening.ScreeningType:
+        """
+        Get the current screening model of the engine.
+        """
+    def getSpeciesDestructionTimescales(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> dict[fourdst._phys.atomic.Species, float]:
+        """
+        Get the destruction timescales for each species in the network.
+        """
+    def getSpeciesIndex(self, species: fourdst._phys.atomic.Species) -> int:
+        """
+        Get the index of a species in the network.
+        """
+    def getSpeciesStatus(self, species: fourdst._phys.atomic.Species) -> SpeciesStatus:
+        """
+        Get the status of a species in the network.
+        """
+    def getSpeciesTimescales(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> dict[fourdst._phys.atomic.Species, float]:
+        """
+        Get the timescales for each species in the network.
+        """
+    def getStoichiometryMatrixEntry(self, species: fourdst._phys.atomic.Species, reaction: ...) -> int:
+        """
+        Get an entry from the stoichiometry matrix.
+        """
+    def isStale(self, netIn: gridfire._gridfire.type.NetIn) -> bool:
+        """
+        Check if the engine is stale based on the provided NetIn object.
+        """
+    def mapNetInToMolarAbundanceVector(self, netIn: gridfire._gridfire.type.NetIn) -> list[float]:
+        """
+        Map a NetIn object to a vector of molar abundances.
+        """
+    def primeEngine(self, netIn: gridfire._gridfire.type.NetIn) -> PrimingReport:
+        """
+        Prime the engine with a NetIn object to prepare for calculations.
+        """
+    def rebuild(self, composition: ..., depth: gridfire._gridfire.engine.NetworkBuildDepth | typing.SupportsInt = ...) -> None:
+        """
+        Rebuild the engine with a new composition and build depth.
+        """
+    def setNetworkReactions(self, reactions: gridfire._gridfire.reaction.ReactionSet) -> None:
+        """
+        Set the network reactions to a new set of reactions.
+        """
+    def setScreeningModel(self, screeningModel: gridfire._gridfire.screening.ScreeningType) -> None:
+        """
+        Set the screening model for the engine.
+        """
+    def update(self, netIn: gridfire._gridfire.type.NetIn) -> fourdst._phys.composition.Composition:
+        """
+        Update the engine state based on the provided NetIn object.
+        """
+class BuildDepthType:
+    pass
+class DefinedEngineView(DynamicEngine):
+    def __init__(self, peNames: collections.abc.Sequence[str], baseEngine: GraphEngine) -> None:
+        """
+        Construct a defined engine view with a list of tracked reactions and a base engine.
+        """
+    def calculateEpsDerivatives(self, comp: ..., T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> ...:
+        """
+        Calculate deps/dT and deps/drho
+        """
+    def calculateMolarReactionFlow(self: DynamicEngine, reaction: ..., comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> float:
+        """
+        Calculate the molar reaction flow for a given reaction.
+        """
+    def calculateRHSAndEnergy(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> StepDerivatives:
+        """
+        Calculate the right-hand side (dY/dt) and energy generation rate.
+        """
+    def collectComposition(self, composition: ..., T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> fourdst._phys.composition.Composition:
+        """
+        Recursively collect composition from current engine and any sub engines if they exist.
+        """
+    @typing.overload
+    def generateJacobianMatrix(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> NetworkJacobian:
+        """
+        Generate the Jacobian matrix for the current state.
+        """
+    @typing.overload
+    def generateJacobianMatrix(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat, activeSpecies: collections.abc.Sequence[fourdst._phys.atomic.Species]) -> NetworkJacobian:
+        """
+        Generate the jacobian matrix only for the subset of the matrix representing the active species.
+        """
+    @typing.overload
+    def generateJacobianMatrix(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat, sparsityPattern: collections.abc.Sequence[tuple[typing.SupportsInt, typing.SupportsInt]]) -> NetworkJacobian:
+        """
+        Generate the jacobian matrix for the given sparsity pattern
+        """
+    def generateStoichiometryMatrix(self) -> None:
+        ...
+    def getBaseEngine(self) -> DynamicEngine:
+        """
+        Get the base engine associated with this defined engine view.
+        """
+    def getDepth(self) -> gridfire._gridfire.engine.NetworkBuildDepth | int:
+        """
+        Get the current build depth of the engine.
+        """
+    def getNetworkReactions(self) -> gridfire._gridfire.reaction.ReactionSet:
+        """
+        Get the set of logical reactions in the network.
+        """
+    def getNetworkSpecies(self) -> list[fourdst._phys.atomic.Species]:
+        """
+        Get the list of species in the network.
+        """
+    def getScreeningModel(self) -> gridfire._gridfire.screening.ScreeningType:
+        """
+        Get the current screening model of the engine.
+        """
+    def getSpeciesDestructionTimescales(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> dict[fourdst._phys.atomic.Species, float]:
+        """
+        Get the destruction timescales for each species in the network.
+        """
+    def getSpeciesIndex(self, species: fourdst._phys.atomic.Species) -> int:
+        """
+        Get the index of a species in the network.
+        """
+    def getSpeciesStatus(self, species: fourdst._phys.atomic.Species) -> SpeciesStatus:
+        """
+        Get the status of a species in the network.
+        """
+    def getSpeciesTimescales(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> dict[fourdst._phys.atomic.Species, float]:
+        """
+        Get the timescales for each species in the network.
+        """
+    def getStoichiometryMatrixEntry(self, species: fourdst._phys.atomic.Species, reaction: ...) -> int:
+        """
+        Get an entry from the stoichiometry matrix.
+        """
+    def isStale(self, netIn: gridfire._gridfire.type.NetIn) -> bool:
+        """
+        Check if the engine is stale based on the provided NetIn object.
+        """
+    def mapNetInToMolarAbundanceVector(self, netIn: gridfire._gridfire.type.NetIn) -> list[float]:
+        """
+        Map a NetIn object to a vector of molar abundances.
+        """
+    def primeEngine(self, netIn: gridfire._gridfire.type.NetIn) -> PrimingReport:
+        """
+        Prime the engine with a NetIn object to prepare for calculations.
+        """
+    def rebuild(self, composition: ..., depth: gridfire._gridfire.engine.NetworkBuildDepth | typing.SupportsInt = ...) -> None:
+        """
+        Rebuild the engine with a new composition and build depth.
+        """
+    def setNetworkReactions(self, reactions: gridfire._gridfire.reaction.ReactionSet) -> None:
+        """
+        Set the network reactions to a new set of reactions.
+        """
+    def setScreeningModel(self, screeningModel: gridfire._gridfire.screening.ScreeningType) -> None:
+        """
+        Set the screening model for the engine.
+        """
+    def update(self, netIn: gridfire._gridfire.type.NetIn) -> fourdst._phys.composition.Composition:
+        """
+        Update the engine state based on the provided NetIn object.
+        """
+class DynamicEngine:
+    pass
+class Engine:
+    pass
+class EngineTypes:
+    """
+    Members:
+    
+      GRAPH_ENGINE : The standard graph-based engine.
+    
+      ADAPTIVE_ENGINE_VIEW : An engine that adapts based on certain criteria.
+    
+      MULTISCALE_PARTITIONING_ENGINE_VIEW : An engine that partitions the system at multiple scales.
+    
+      PRIMING_ENGINE_VIEW : An engine that uses a priming strategy for simulations.
+    
+      DEFINED_ENGINE_VIEW : An engine defined by user specifications.
+    
+      FILE_DEFINED_ENGINE_VIEW : An engine defined through external files.
+    """
+    ADAPTIVE_ENGINE_VIEW: typing.ClassVar[EngineTypes]  # value = <EngineTypes.ADAPTIVE_ENGINE_VIEW: 1>
+    DEFINED_ENGINE_VIEW: typing.ClassVar[EngineTypes]  # value = <EngineTypes.DEFINED_ENGINE_VIEW: 4>
+    FILE_DEFINED_ENGINE_VIEW: typing.ClassVar[EngineTypes]  # value = <EngineTypes.FILE_DEFINED_ENGINE_VIEW: 5>
+    GRAPH_ENGINE: typing.ClassVar[EngineTypes]  # value = <EngineTypes.GRAPH_ENGINE: 0>
+    MULTISCALE_PARTITIONING_ENGINE_VIEW: typing.ClassVar[EngineTypes]  # value = <EngineTypes.MULTISCALE_PARTITIONING_ENGINE_VIEW: 2>
+    PRIMING_ENGINE_VIEW: typing.ClassVar[EngineTypes]  # value = <EngineTypes.PRIMING_ENGINE_VIEW: 3>
+    __members__: typing.ClassVar[dict[str, EngineTypes]]  # value = {'GRAPH_ENGINE': <EngineTypes.GRAPH_ENGINE: 0>, 'ADAPTIVE_ENGINE_VIEW': <EngineTypes.ADAPTIVE_ENGINE_VIEW: 1>, 'MULTISCALE_PARTITIONING_ENGINE_VIEW': <EngineTypes.MULTISCALE_PARTITIONING_ENGINE_VIEW: 2>, 'PRIMING_ENGINE_VIEW': <EngineTypes.PRIMING_ENGINE_VIEW: 3>, 'DEFINED_ENGINE_VIEW': <EngineTypes.DEFINED_ENGINE_VIEW: 4>, 'FILE_DEFINED_ENGINE_VIEW': <EngineTypes.FILE_DEFINED_ENGINE_VIEW: 5>}
+    def __eq__(self, other: typing.Any) -> bool:
+        ...
+    def __getstate__(self) -> int:
+        ...
+    def __hash__(self) -> int:
+        ...
+    def __index__(self) -> int:
+        ...
+    def __init__(self, value: typing.SupportsInt) -> None:
+        ...
+    def __int__(self) -> int:
+        ...
+    def __ne__(self, other: typing.Any) -> bool:
+        ...
+    @typing.overload
+    def __repr__(self) -> str:
+        ...
+    @typing.overload
+    def __repr__(self) -> str:
+        """
+        String representation of the EngineTypes.
+        """
+    def __setstate__(self, state: typing.SupportsInt) -> None:
+        ...
+    def __str__(self) -> str:
+        ...
+    @property
+    def name(self) -> str:
+        ...
+    @property
+    def value(self) -> int:
+        ...
+class FileDefinedEngineView(DefinedEngineView):
+    def __init__(self, baseEngine: GraphEngine, fileName: str, parser: gridfire._gridfire.io.NetworkFileParser) -> None:
+        """
+        Construct a defined engine view from a file and a base engine.
+        """
+    def calculateEpsDerivatives(self, comp: ..., T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> ...:
+        """
+        Calculate deps/dT and deps/drho
+        """
+    def calculateMolarReactionFlow(self: DynamicEngine, reaction: ..., comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> float:
+        """
+        Calculate the molar reaction flow for a given reaction.
+        """
+    def calculateRHSAndEnergy(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> StepDerivatives:
+        """
+        Calculate the right-hand side (dY/dt) and energy generation rate.
+        """
+    def collectComposition(self, composition: ..., T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> fourdst._phys.composition.Composition:
+        """
+        Recursively collect composition from current engine and any sub engines if they exist.
+        """
+    @typing.overload
+    def generateJacobianMatrix(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> NetworkJacobian:
+        """
+        Generate the Jacobian matrix for the current state.
+        """
+    @typing.overload
+    def generateJacobianMatrix(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat, activeSpecies: collections.abc.Sequence[fourdst._phys.atomic.Species]) -> NetworkJacobian:
+        """
+        Generate the jacobian matrix only for the subset of the matrix representing the active species.
+        """
+    @typing.overload
+    def generateJacobianMatrix(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat, sparsityPattern: collections.abc.Sequence[tuple[typing.SupportsInt, typing.SupportsInt]]) -> NetworkJacobian:
+        """
+        Generate the jacobian matrix for the given sparsity pattern
+        """
+    def generateStoichiometryMatrix(self) -> None:
+        ...
+    def getBaseEngine(self) -> DynamicEngine:
+        """
+        Get the base engine associated with this file defined engine view.
+        """
+    def getDepth(self) -> gridfire._gridfire.engine.NetworkBuildDepth | int:
+        """
+        Get the current build depth of the engine.
+        """
+    def getNetworkFile(self) -> str:
+        """
+        Get the network file associated with this defined engine view.
+        """
+    def getNetworkReactions(self) -> gridfire._gridfire.reaction.ReactionSet:
+        """
+        Get the set of logical reactions in the network.
+        """
+    def getNetworkSpecies(self) -> list[fourdst._phys.atomic.Species]:
+        """
+        Get the list of species in the network.
+        """
+    def getParser(self) -> gridfire._gridfire.io.NetworkFileParser:
+        """
+        Get the parser used for this defined engine view.
+        """
+    def getScreeningModel(self) -> gridfire._gridfire.screening.ScreeningType:
+        """
+        Get the current screening model of the engine.
+        """
+    def getSpeciesDestructionTimescales(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> dict[fourdst._phys.atomic.Species, float]:
+        """
+        Get the destruction timescales for each species in the network.
+        """
+    def getSpeciesIndex(self, species: fourdst._phys.atomic.Species) -> int:
+        """
+        Get the index of a species in the network.
+        """
+    def getSpeciesStatus(self, species: fourdst._phys.atomic.Species) -> SpeciesStatus:
+        """
+        Get the status of a species in the network.
+        """
+    def getSpeciesTimescales(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> dict[fourdst._phys.atomic.Species, float]:
+        """
+        Get the timescales for each species in the network.
+        """
+    def getStoichiometryMatrixEntry(self, species: fourdst._phys.atomic.Species, reaction: ...) -> int:
+        """
+        Get an entry from the stoichiometry matrix.
+        """
+    def isStale(self, netIn: gridfire._gridfire.type.NetIn) -> bool:
+        """
+        Check if the engine is stale based on the provided NetIn object.
+        """
+    def mapNetInToMolarAbundanceVector(self, netIn: gridfire._gridfire.type.NetIn) -> list[float]:
+        """
+        Map a NetIn object to a vector of molar abundances.
+        """
+    def primeEngine(self, netIn: gridfire._gridfire.type.NetIn) -> PrimingReport:
+        """
+        Prime the engine with a NetIn object to prepare for calculations.
+        """
+    def rebuild(self, composition: ..., depth: gridfire._gridfire.engine.NetworkBuildDepth | typing.SupportsInt = ...) -> None:
+        """
+        Rebuild the engine with a new composition and build depth.
+        """
+    def setNetworkReactions(self, reactions: gridfire._gridfire.reaction.ReactionSet) -> None:
+        """
+        Set the network reactions to a new set of reactions.
+        """
+    def setScreeningModel(self, screeningModel: gridfire._gridfire.screening.ScreeningType) -> None:
+        """
+        Set the screening model for the engine.
+        """
+    def update(self, netIn: gridfire._gridfire.type.NetIn) -> fourdst._phys.composition.Composition:
+        """
+        Update the engine state based on the provided NetIn object.
+        """
+class GraphEngine(DynamicEngine):
+    @staticmethod
+    def getNetReactionStoichiometry(reaction: ...) -> dict[fourdst._phys.atomic.Species, int]:
+        """
+        Get the net stoichiometry for a given reaction.
+        """
+    @typing.overload
+    def __init__(self, composition: fourdst._phys.composition.Composition, depth: gridfire._gridfire.engine.NetworkBuildDepth | typing.SupportsInt = ...) -> None:
+        """
+        Initialize GraphEngine with a composition and build depth.
+        """
+    @typing.overload
+    def __init__(self, composition: fourdst._phys.composition.Composition, partitionFunction: gridfire._gridfire.partition.PartitionFunction, depth: gridfire._gridfire.engine.NetworkBuildDepth | typing.SupportsInt = ...) -> None:
+        """
+        Initialize GraphEngine with a composition, partition function and build depth.
+        """
+    @typing.overload
+    def __init__(self, reactions: gridfire._gridfire.reaction.ReactionSet) -> None:
+        """
+        Initialize GraphEngine with a set of reactions.
+        """
+    def calculateEpsDerivatives(self, comp: ..., T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> ...:
+        """
+        Calculate deps/dT and deps/drho
+        """
+    def calculateMolarReactionFlow(self: DynamicEngine, reaction: ..., comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> float:
+        """
+        Calculate the molar reaction flow for a given reaction.
+        """
+    def calculateRHSAndEnergy(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> StepDerivatives:
+        """
+        Calculate the right-hand side (dY/dt) and energy generation rate.
+        """
+    def calculateReverseRate(self, reaction: ..., T9: typing.SupportsFloat, rho: typing.SupportsFloat, composition: ...) -> float:
+        """
+        Calculate the reverse rate for a given reaction at a specific temperature, density, and composition.
+        """
+    def calculateReverseRateTwoBody(self, reaction: ..., T9: typing.SupportsFloat, forwardRate: typing.SupportsFloat, expFactor: typing.SupportsFloat) -> float:
+        """
+        Calculate the reverse rate for a two-body reaction at a specific temperature.
+        """
+    def calculateReverseRateTwoBodyDerivative(self, reaction: ..., T9: typing.SupportsFloat, rho: typing.SupportsFloat, composition: fourdst._phys.composition.Composition, reverseRate: typing.SupportsFloat) -> float:
+        """
+        Calculate the derivative of the reverse rate for a two-body reaction at a specific temperature.
+        """
+    def collectComposition(self, composition: ..., T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> fourdst._phys.composition.Composition:
+        """
+        Recursively collect composition from current engine and any sub engines if they exist.
+        """
+    def exportToCSV(self, filename: str) -> None:
+        """
+        Export the network to a CSV file for analysis.
+        """
+    def exportToDot(self, filename: str) -> None:
+        """
+        Export the network to a DOT file for visualization.
+        """
+    @typing.overload
+    def generateJacobianMatrix(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> NetworkJacobian:
+        """
+        Generate the Jacobian matrix for the current state.
+        """
+    @typing.overload
+    def generateJacobianMatrix(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat, activeSpecies: collections.abc.Sequence[fourdst._phys.atomic.Species]) -> NetworkJacobian:
+        """
+        Generate the jacobian matrix only for the subset of the matrix representing the active species.
+        """
+    @typing.overload
+    def generateJacobianMatrix(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat, sparsityPattern: collections.abc.Sequence[tuple[typing.SupportsInt, typing.SupportsInt]]) -> NetworkJacobian:
+        """
+        Generate the jacobian matrix for the given sparsity pattern
+        """
+    def generateStoichiometryMatrix(self) -> None:
+        ...
+    def getDepth(self) -> gridfire._gridfire.engine.NetworkBuildDepth | int:
+        """
+        Get the current build depth of the engine.
+        """
+    def getNetworkReactions(self) -> gridfire._gridfire.reaction.ReactionSet:
+        """
+        Get the set of logical reactions in the network.
+        """
+    def getNetworkSpecies(self) -> list[fourdst._phys.atomic.Species]:
+        """
+        Get the list of species in the network.
+        """
+    def getPartitionFunction(self) -> gridfire._gridfire.partition.PartitionFunction:
+        """
+        Get the partition function used by the engine.
+        """
+    def getScreeningModel(self) -> gridfire._gridfire.screening.ScreeningType:
+        """
+        Get the current screening model of the engine.
+        """
+    @typing.overload
+    def getSpeciesDestructionTimescales(self, composition: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat, activeReactions: gridfire._gridfire.reaction.ReactionSet) -> ...:
+        ...
+    @typing.overload
+    def getSpeciesDestructionTimescales(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> dict[fourdst._phys.atomic.Species, float]:
+        """
+        Get the destruction timescales for each species in the network.
+        """
+    def getSpeciesIndex(self, species: fourdst._phys.atomic.Species) -> int:
+        """
+        Get the index of a species in the network.
+        """
+    def getSpeciesStatus(self, species: fourdst._phys.atomic.Species) -> SpeciesStatus:
+        """
+        Get the status of a species in the network.
+        """
+    @typing.overload
+    def getSpeciesTimescales(self, composition: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat, activeReactions: gridfire._gridfire.reaction.ReactionSet) -> ...:
+        ...
+    @typing.overload
+    def getSpeciesTimescales(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> dict[fourdst._phys.atomic.Species, float]:
+        """
+        Get the timescales for each species in the network.
+        """
+    def getStoichiometryMatrixEntry(self, species: fourdst._phys.atomic.Species, reaction: ...) -> int:
+        """
+        Get an entry from the stoichiometry matrix.
+        """
+    def involvesSpecies(self, species: fourdst._phys.atomic.Species) -> bool:
+        """
+        Check if a given species is involved in the network.
+        """
+    def isPrecomputationEnabled(self) -> bool:
+        """
+        Check if precomputation is enabled for the engine.
+        """
+    def isStale(self, netIn: gridfire._gridfire.type.NetIn) -> bool:
+        """
+        Check if the engine is stale based on the provided NetIn object.
+        """
+    def isUsingReverseReactions(self) -> bool:
+        """
+        Check if the engine is using reverse reactions.
+        """
+    def mapNetInToMolarAbundanceVector(self, netIn: gridfire._gridfire.type.NetIn) -> list[float]:
+        """
+        Map a NetIn object to a vector of molar abundances.
+        """
+    def primeEngine(self, netIn: gridfire._gridfire.type.NetIn) -> PrimingReport:
+        """
+        Prime the engine with a NetIn object to prepare for calculations.
+        """
+    def rebuild(self, composition: ..., depth: gridfire._gridfire.engine.NetworkBuildDepth | typing.SupportsInt = ...) -> None:
+        """
+        Rebuild the engine with a new composition and build depth.
+        """
+    def setNetworkReactions(self, reactions: gridfire._gridfire.reaction.ReactionSet) -> None:
+        """
+        Set the network reactions to a new set of reactions.
+        """
+    def setPrecomputation(self, precompute: bool) -> None:
+        """
+        Enable or disable precomputation for the engine.
+        """
+    def setScreeningModel(self, screeningModel: gridfire._gridfire.screening.ScreeningType) -> None:
+        """
+        Set the screening model for the engine.
+        """
+    def setUseReverseReactions(self, useReverse: bool) -> None:
+        """
+        Enable or disable the use of reverse reactions in the engine.
+        """
+    def update(self, netIn: gridfire._gridfire.type.NetIn) -> fourdst._phys.composition.Composition:
+        """
+        Update the engine state based on the provided NetIn object.
+        """
+class MultiscalePartitioningEngineView(DynamicEngine):
+    def __init__(self, baseEngine: GraphEngine) -> None:
+        """
+        Construct a multiscale partitioning engine view with a base engine.
+        """
+    def calculateEpsDerivatives(self, comp: ..., T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> ...:
+        """
+        Calculate deps/dT and deps/drho
+        """
+    def calculateMolarReactionFlow(self: DynamicEngine, reaction: ..., comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> float:
+        """
+        Calculate the molar reaction flow for a given reaction.
+        """
+    def calculateRHSAndEnergy(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> StepDerivatives:
+        """
+        Calculate the right-hand side (dY/dt) and energy generation rate.
+        """
+    def collectComposition(self, composition: ..., T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> fourdst._phys.composition.Composition:
+        """
+        Recursively collect composition from current engine and any sub engines if they exist.
+        """
+    def exportToDot(self, filename: str, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> None:
+        """
+        Export the network to a DOT file for visualization.
+        """
+    @typing.overload
+    def generateJacobianMatrix(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> NetworkJacobian:
+        """
+        Generate the Jacobian matrix for the current state.
+        """
+    @typing.overload
+    def generateJacobianMatrix(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat, activeSpecies: collections.abc.Sequence[fourdst._phys.atomic.Species]) -> NetworkJacobian:
+        """
+        Generate the jacobian matrix only for the subset of the matrix representing the active species.
+        """
+    @typing.overload
+    def generateJacobianMatrix(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat, sparsityPattern: collections.abc.Sequence[tuple[typing.SupportsInt, typing.SupportsInt]]) -> NetworkJacobian:
+        """
+        Generate the jacobian matrix for the given sparsity pattern
+        """
+    def generateStoichiometryMatrix(self) -> None:
+        ...
+    def getBaseEngine(self) -> DynamicEngine:
+        """
+        Get the base engine associated with this multiscale partitioning engine view.
+        """
+    def getDepth(self) -> gridfire._gridfire.engine.NetworkBuildDepth | int:
+        """
+        Get the current build depth of the engine.
+        """
+    def getDynamicSpecies(self) -> list[fourdst._phys.atomic.Species]:
+        """
+        Get the list of dynamic species in the network.
+        """
+    def getFastSpecies(self) -> list[fourdst._phys.atomic.Species]:
+        """
+        Get the list of fast species in the network.
+        """
+    def getNetworkReactions(self) -> gridfire._gridfire.reaction.ReactionSet:
+        """
+        Get the set of logical reactions in the network.
+        """
+    def getNetworkSpecies(self) -> list[fourdst._phys.atomic.Species]:
+        """
+        Get the list of species in the network.
+        """
+    def getNormalizedEquilibratedComposition(self, comp: ..., T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> fourdst._phys.composition.Composition:
+        """
+        Get the normalized equilibrated composition for the algebraic species.
+        """
+    def getScreeningModel(self) -> gridfire._gridfire.screening.ScreeningType:
+        """
+        Get the current screening model of the engine.
+        """
+    def getSpeciesDestructionTimescales(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> dict[fourdst._phys.atomic.Species, float]:
+        """
+        Get the destruction timescales for each species in the network.
+        """
+    def getSpeciesIndex(self, species: fourdst._phys.atomic.Species) -> int:
+        """
+        Get the index of a species in the network.
+        """
+    def getSpeciesStatus(self, species: fourdst._phys.atomic.Species) -> SpeciesStatus:
+        """
+        Get the status of a species in the network.
+        """
+    def getSpeciesTimescales(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> dict[fourdst._phys.atomic.Species, float]:
+        """
+        Get the timescales for each species in the network.
+        """
+    def getStoichiometryMatrixEntry(self, species: fourdst._phys.atomic.Species, reaction: ...) -> int:
+        """
+        Get an entry from the stoichiometry matrix.
+        """
+    def involvesSpecies(self, species: fourdst._phys.atomic.Species) -> bool:
+        """
+        Check if a given species is involved in the network (in either the algebraic or dynamic set).
+        """
+    def involvesSpeciesInDynamic(self, species: fourdst._phys.atomic.Species) -> bool:
+        """
+        Check if a given species is involved in the network's dynamic set.
+        """
+    def involvesSpeciesInQSE(self, species: fourdst._phys.atomic.Species) -> bool:
+        """
+        Check if a given species is involved in the network's algebraic set.
+        """
+    def isStale(self, netIn: gridfire._gridfire.type.NetIn) -> bool:
+        """
+        Check if the engine is stale based on the provided NetIn object.
+        """
+    def mapNetInToMolarAbundanceVector(self, netIn: gridfire._gridfire.type.NetIn) -> list[float]:
+        """
+        Map a NetIn object to a vector of molar abundances.
+        """
+    @typing.overload
+    def partitionNetwork(self, netIn: gridfire._gridfire.type.NetIn) -> fourdst._phys.composition.Composition:
+        """
+        Partition the network based on species timescales and connectivity.
+        """
+    @typing.overload
+    def partitionNetwork(self, netIn: gridfire._gridfire.type.NetIn) -> fourdst._phys.composition.Composition:
+        """
+        Partition the network based on a NetIn object.
+        """
+    def primeEngine(self, netIn: gridfire._gridfire.type.NetIn) -> PrimingReport:
+        """
+        Prime the engine with a NetIn object to prepare for calculations.
+        """
+    def rebuild(self, composition: ..., depth: gridfire._gridfire.engine.NetworkBuildDepth | typing.SupportsInt = ...) -> None:
+        """
+        Rebuild the engine with a new composition and build depth.
+        """
+    def setNetworkReactions(self, reactions: gridfire._gridfire.reaction.ReactionSet) -> None:
+        """
+        Set the network reactions to a new set of reactions.
+        """
+    def setScreeningModel(self, screeningModel: gridfire._gridfire.screening.ScreeningType) -> None:
+        """
+        Set the screening model for the engine.
+        """
+    def update(self, netIn: gridfire._gridfire.type.NetIn) -> fourdst._phys.composition.Composition:
+        """
+        Update the engine state based on the provided NetIn object.
+        """
+class NetworkBuildDepth:
+    """
+    Members:
+    
+      Full : Full network build depth
+    
+      Shallow : Shallow network build depth
+    
+      SecondOrder : Second order network build depth
+    
+      ThirdOrder : Third order network build depth
+    
+      FourthOrder : Fourth order network build depth
+    
+      FifthOrder : Fifth order network build depth
+    """
+    FifthOrder: typing.ClassVar[NetworkBuildDepth]  # value = <NetworkBuildDepth.FifthOrder: 5>
+    FourthOrder: typing.ClassVar[NetworkBuildDepth]  # value = <NetworkBuildDepth.FourthOrder: 4>
+    Full: typing.ClassVar[NetworkBuildDepth]  # value = <NetworkBuildDepth.Full: -1>
+    SecondOrder: typing.ClassVar[NetworkBuildDepth]  # value = <NetworkBuildDepth.SecondOrder: 2>
+    Shallow: typing.ClassVar[NetworkBuildDepth]  # value = <NetworkBuildDepth.Shallow: 1>
+    ThirdOrder: typing.ClassVar[NetworkBuildDepth]  # value = <NetworkBuildDepth.ThirdOrder: 3>
+    __members__: typing.ClassVar[dict[str, NetworkBuildDepth]]  # value = {'Full': <NetworkBuildDepth.Full: -1>, 'Shallow': <NetworkBuildDepth.Shallow: 1>, 'SecondOrder': <NetworkBuildDepth.SecondOrder: 2>, 'ThirdOrder': <NetworkBuildDepth.ThirdOrder: 3>, 'FourthOrder': <NetworkBuildDepth.FourthOrder: 4>, 'FifthOrder': <NetworkBuildDepth.FifthOrder: 5>}
+    def __eq__(self, other: typing.Any) -> bool:
+        ...
+    def __getstate__(self) -> int:
+        ...
+    def __hash__(self) -> int:
+        ...
+    def __index__(self) -> int:
+        ...
+    def __init__(self, value: typing.SupportsInt) -> None:
+        ...
+    def __int__(self) -> int:
+        ...
+    def __ne__(self, other: typing.Any) -> bool:
+        ...
+    def __repr__(self) -> str:
+        ...
+    def __setstate__(self, state: typing.SupportsInt) -> None:
+        ...
+    def __str__(self) -> str:
+        ...
+    @property
+    def name(self) -> str:
+        ...
+    @property
+    def value(self) -> int:
+        ...
+class NetworkConstructionFlags:
+    """
+    Members:
+    
+      NONE : No special construction flags.
+    
+      REACLIB_STRONG : Include strong reactions from reaclib.
+    
+      WRL_BETA_MINUS : Include beta-minus decay reactions from weak rate library.
+    
+      WRL_BETA_PLUS : Include beta-plus decay reactions from weak rate library.
+    
+      WRL_ELECTRON_CAPTURE : Include electron capture reactions from weak rate library.
+    
+      WRL_POSITRON_CAPTURE : Include positron capture reactions from weak rate library.
+    
+      REACLIB_WEAK : Include weak reactions from reaclib.
+    
+      WRL_WEAK : Include all weak reactions from weak rate library.
+    
+      REACLIB : Include all reactions from reaclib.
+    
+      DEFAULT : Default construction flags (Reaclib strong and weak).
+    """
+    DEFAULT: typing.ClassVar[NetworkConstructionFlags]  # value = <NetworkConstructionFlags.REACLIB: 33>
+    NONE: typing.ClassVar[NetworkConstructionFlags]  # value = <NetworkConstructionFlags.NONE: 0>
+    REACLIB: typing.ClassVar[NetworkConstructionFlags]  # value = <NetworkConstructionFlags.REACLIB: 33>
+    REACLIB_STRONG: typing.ClassVar[NetworkConstructionFlags]  # value = <NetworkConstructionFlags.REACLIB_STRONG: 1>
+    REACLIB_WEAK: typing.ClassVar[NetworkConstructionFlags]  # value = <NetworkConstructionFlags.REACLIB_WEAK: 32>
+    WRL_BETA_MINUS: typing.ClassVar[NetworkConstructionFlags]  # value = <NetworkConstructionFlags.WRL_BETA_MINUS: 2>
+    WRL_BETA_PLUS: typing.ClassVar[NetworkConstructionFlags]  # value = <NetworkConstructionFlags.WRL_BETA_PLUS: 4>
+    WRL_ELECTRON_CAPTURE: typing.ClassVar[NetworkConstructionFlags]  # value = <NetworkConstructionFlags.WRL_ELECTRON_CAPTURE: 8>
+    WRL_POSITRON_CAPTURE: typing.ClassVar[NetworkConstructionFlags]  # value = <NetworkConstructionFlags.WRL_POSITRON_CAPTURE: 16>
+    WRL_WEAK: typing.ClassVar[NetworkConstructionFlags]  # value = <NetworkConstructionFlags.WRL_WEAK: 30>
+    __members__: typing.ClassVar[dict[str, NetworkConstructionFlags]]  # value = {'NONE': <NetworkConstructionFlags.NONE: 0>, 'REACLIB_STRONG': <NetworkConstructionFlags.REACLIB_STRONG: 1>, 'WRL_BETA_MINUS': <NetworkConstructionFlags.WRL_BETA_MINUS: 2>, 'WRL_BETA_PLUS': <NetworkConstructionFlags.WRL_BETA_PLUS: 4>, 'WRL_ELECTRON_CAPTURE': <NetworkConstructionFlags.WRL_ELECTRON_CAPTURE: 8>, 'WRL_POSITRON_CAPTURE': <NetworkConstructionFlags.WRL_POSITRON_CAPTURE: 16>, 'REACLIB_WEAK': <NetworkConstructionFlags.REACLIB_WEAK: 32>, 'WRL_WEAK': <NetworkConstructionFlags.WRL_WEAK: 30>, 'REACLIB': <NetworkConstructionFlags.REACLIB: 33>, 'DEFAULT': <NetworkConstructionFlags.REACLIB: 33>}
+    def __eq__(self, other: typing.Any) -> bool:
+        ...
+    def __getstate__(self) -> int:
+        ...
+    def __hash__(self) -> int:
+        ...
+    def __index__(self) -> int:
+        ...
+    def __init__(self, value: typing.SupportsInt) -> None:
+        ...
+    def __int__(self) -> int:
+        ...
+    def __ne__(self, other: typing.Any) -> bool:
+        ...
+    @typing.overload
+    def __repr__(self) -> str:
+        ...
+    @typing.overload
+    def __repr__(self) -> str:
+        ...
+    def __setstate__(self, state: typing.SupportsInt) -> None:
+        ...
+    def __str__(self) -> str:
+        ...
+    @property
+    def name(self) -> str:
+        ...
+    @property
+    def value(self) -> int:
+        ...
+class NetworkJacobian:
+    @typing.overload
+    def __getitem__(self, key: tuple[fourdst._phys.atomic.Species, fourdst._phys.atomic.Species]) -> float:
+        """
+        Get an entry from the Jacobian matrix using species identifiers.
+        """
+    @typing.overload
+    def __getitem__(self, key: tuple[typing.SupportsInt, typing.SupportsInt]) -> float:
+        """
+        Get an entry from the Jacobian matrix using indices.
+        """
+    @typing.overload
+    def __setitem__(self, key: tuple[fourdst._phys.atomic.Species, fourdst._phys.atomic.Species], value: typing.SupportsFloat) -> None:
+        """
+        Set an entry in the Jacobian matrix using species identifiers.
+        """
+    @typing.overload
+    def __setitem__(self, key: tuple[typing.SupportsInt, typing.SupportsInt], value: typing.SupportsFloat) -> None:
+        """
+        Set an entry in the Jacobian matrix using indices.
+        """
+    def data(self) -> ...:
+        """
+        Get the underlying sparse matrix data.
+        """
+    def infs(self) -> list[tuple[tuple[fourdst._phys.atomic.Species, fourdst._phys.atomic.Species], float]]:
+        """
+        Get all infinite entries in the Jacobian matrix.
+        """
+    def mapping(self) -> dict[fourdst._phys.atomic.Species, int]:
+        """
+        Get the species-to-index mapping.
+        """
+    def nans(self) -> list[tuple[tuple[fourdst._phys.atomic.Species, fourdst._phys.atomic.Species], float]]:
+        """
+        Get all NaN entries in the Jacobian matrix.
+        """
+    def nnz(self) -> int:
+        """
+        Get the number of non-zero entries in the Jacobian matrix.
+        """
+    def rank(self) -> int:
+        """
+        Get the rank of the Jacobian matrix.
+        """
+    def shape(self) -> tuple[int, int]:
+        """
+        Get the shape of the Jacobian matrix as (rows, columns).
+        """
+    def singular(self) -> bool:
+        """
+        Check if the Jacobian matrix is singular.
+        """
+    def to_csv(self, filename: str) -> None:
+        """
+        Export the Jacobian matrix to a CSV file.
+        """
+    def to_numpy(self) -> numpy.typing.NDArray[numpy.float64]:
+        """
+        Convert the Jacobian matrix to a NumPy array.
+        """
+class NetworkPrimingEngineView(DefinedEngineView):
+    @typing.overload
+    def __init__(self, primingSymbol: str, baseEngine: GraphEngine) -> None:
+        """
+        Construct a priming engine view with a priming symbol and a base engine.
+        """
+    @typing.overload
+    def __init__(self, primingSpecies: fourdst._phys.atomic.Species, baseEngine: GraphEngine) -> None:
+        """
+        Construct a priming engine view with a priming species and a base engine.
+        """
+    def calculateEpsDerivatives(self, comp: ..., T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> ...:
+        """
+        Calculate deps/dT and deps/drho
+        """
+    def calculateMolarReactionFlow(self: DynamicEngine, reaction: ..., comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> float:
+        """
+        Calculate the molar reaction flow for a given reaction.
+        """
+    def calculateRHSAndEnergy(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> StepDerivatives:
+        """
+        Calculate the right-hand side (dY/dt) and energy generation rate.
+        """
+    def collectComposition(self, composition: ..., T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> fourdst._phys.composition.Composition:
+        """
+        Recursively collect composition from current engine and any sub engines if they exist.
+        """
+    @typing.overload
+    def generateJacobianMatrix(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> NetworkJacobian:
+        """
+        Generate the Jacobian matrix for the current state.
+        """
+    @typing.overload
+    def generateJacobianMatrix(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat, activeSpecies: collections.abc.Sequence[fourdst._phys.atomic.Species]) -> NetworkJacobian:
+        """
+        Generate the jacobian matrix only for the subset of the matrix representing the active species.
+        """
+    @typing.overload
+    def generateJacobianMatrix(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat, sparsityPattern: collections.abc.Sequence[tuple[typing.SupportsInt, typing.SupportsInt]]) -> NetworkJacobian:
+        """
+        Generate the jacobian matrix for the given sparsity pattern
+        """
+    def generateStoichiometryMatrix(self) -> None:
+        ...
+    def getBaseEngine(self) -> DynamicEngine:
+        """
+        Get the base engine associated with this priming engine view.
+        """
+    def getDepth(self) -> gridfire._gridfire.engine.NetworkBuildDepth | int:
+        """
+        Get the current build depth of the engine.
+        """
+    def getNetworkReactions(self) -> gridfire._gridfire.reaction.ReactionSet:
+        """
+        Get the set of logical reactions in the network.
+        """
+    def getNetworkSpecies(self) -> list[fourdst._phys.atomic.Species]:
+        """
+        Get the list of species in the network.
+        """
+    def getScreeningModel(self) -> gridfire._gridfire.screening.ScreeningType:
+        """
+        Get the current screening model of the engine.
+        """
+    def getSpeciesDestructionTimescales(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> dict[fourdst._phys.atomic.Species, float]:
+        """
+        Get the destruction timescales for each species in the network.
+        """
+    def getSpeciesIndex(self, species: fourdst._phys.atomic.Species) -> int:
+        """
+        Get the index of a species in the network.
+        """
+    def getSpeciesStatus(self, species: fourdst._phys.atomic.Species) -> SpeciesStatus:
+        """
+        Get the status of a species in the network.
+        """
+    def getSpeciesTimescales(self: DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> dict[fourdst._phys.atomic.Species, float]:
+        """
+        Get the timescales for each species in the network.
+        """
+    def getStoichiometryMatrixEntry(self, species: fourdst._phys.atomic.Species, reaction: ...) -> int:
+        """
+        Get an entry from the stoichiometry matrix.
+        """
+    def isStale(self, netIn: gridfire._gridfire.type.NetIn) -> bool:
+        """
+        Check if the engine is stale based on the provided NetIn object.
+        """
+    def mapNetInToMolarAbundanceVector(self, netIn: gridfire._gridfire.type.NetIn) -> list[float]:
+        """
+        Map a NetIn object to a vector of molar abundances.
+        """
+    def primeEngine(self, netIn: gridfire._gridfire.type.NetIn) -> PrimingReport:
+        """
+        Prime the engine with a NetIn object to prepare for calculations.
+        """
+    def rebuild(self, composition: ..., depth: gridfire._gridfire.engine.NetworkBuildDepth | typing.SupportsInt = ...) -> None:
+        """
+        Rebuild the engine with a new composition and build depth.
+        """
+    def setNetworkReactions(self, reactions: gridfire._gridfire.reaction.ReactionSet) -> None:
+        """
+        Set the network reactions to a new set of reactions.
+        """
+    def setScreeningModel(self, screeningModel: gridfire._gridfire.screening.ScreeningType) -> None:
+        """
+        Set the screening model for the engine.
+        """
+    def update(self, netIn: gridfire._gridfire.type.NetIn) -> fourdst._phys.composition.Composition:
+        """
+        Update the engine state based on the provided NetIn object.
+        """
+class PrimingReport:
+    def __repr__(self) -> str:
+        ...
+    @property
+    def primedComposition(self) -> fourdst._phys.composition.Composition:
+        """
+        The composition after priming.
+        """
+    @property
+    def status(self) -> PrimingReportStatus:
+        """
+        Status message from the priming process.
+        """
+    @property
+    def success(self) -> bool:
+        """
+        Indicates if the priming was successful.
+        """
+class PrimingReportStatus:
+    """
+    Members:
+    
+      FULL_SUCCESS : Priming was full successful.
+    
+      NO_SPECIES_TO_PRIME : Solver Failed to converge during priming.
+    
+      MAX_ITERATIONS_REACHED : Engine has already been primed.
+    """
+    FULL_SUCCESS: typing.ClassVar[PrimingReportStatus]  # value = <PrimingReportStatus.FULL_SUCCESS: 0>
+    MAX_ITERATIONS_REACHED: typing.ClassVar[PrimingReportStatus]  # value = <PrimingReportStatus.MAX_ITERATIONS_REACHED: 1>
+    NO_SPECIES_TO_PRIME: typing.ClassVar[PrimingReportStatus]  # value = <PrimingReportStatus.NO_SPECIES_TO_PRIME: 2>
+    __members__: typing.ClassVar[dict[str, PrimingReportStatus]]  # value = {'FULL_SUCCESS': <PrimingReportStatus.FULL_SUCCESS: 0>, 'NO_SPECIES_TO_PRIME': <PrimingReportStatus.NO_SPECIES_TO_PRIME: 2>, 'MAX_ITERATIONS_REACHED': <PrimingReportStatus.MAX_ITERATIONS_REACHED: 1>}
+    def __eq__(self, other: typing.Any) -> bool:
+        ...
+    def __getstate__(self) -> int:
+        ...
+    def __hash__(self) -> int:
+        ...
+    def __index__(self) -> int:
+        ...
+    def __init__(self, value: typing.SupportsInt) -> None:
+        ...
+    def __int__(self) -> int:
+        ...
+    def __ne__(self, other: typing.Any) -> bool:
+        ...
+    @typing.overload
+    def __repr__(self) -> str:
+        ...
+    @typing.overload
+    def __repr__(self) -> str:
+        """
+        String representation of the PrimingReport.
+        """
+    def __setstate__(self, state: typing.SupportsInt) -> None:
+        ...
+    def __str__(self) -> str:
+        ...
+    @property
+    def name(self) -> str:
+        ...
+    @property
+    def value(self) -> int:
+        ...
+class SparsityPattern:
+    pass
+class SpeciesStatus:
+    """
+    Members:
+    
+      ACTIVE : Species is active in the network.
+    
+      EQUILIBRIUM : Species is in equilibrium.
+    
+      INACTIVE_FLOW : Species is inactive due to flow.
+    
+      NOT_PRESENT : Species is not present in the network.
+    """
+    ACTIVE: typing.ClassVar[SpeciesStatus]  # value = <SpeciesStatus.ACTIVE: 0>
+    EQUILIBRIUM: typing.ClassVar[SpeciesStatus]  # value = <SpeciesStatus.EQUILIBRIUM: 1>
+    INACTIVE_FLOW: typing.ClassVar[SpeciesStatus]  # value = <SpeciesStatus.INACTIVE_FLOW: 2>
+    NOT_PRESENT: typing.ClassVar[SpeciesStatus]  # value = <SpeciesStatus.NOT_PRESENT: 3>
+    __members__: typing.ClassVar[dict[str, SpeciesStatus]]  # value = {'ACTIVE': <SpeciesStatus.ACTIVE: 0>, 'EQUILIBRIUM': <SpeciesStatus.EQUILIBRIUM: 1>, 'INACTIVE_FLOW': <SpeciesStatus.INACTIVE_FLOW: 2>, 'NOT_PRESENT': <SpeciesStatus.NOT_PRESENT: 3>}
+    def __eq__(self, other: typing.Any) -> bool:
+        ...
+    def __getstate__(self) -> int:
+        ...
+    def __hash__(self) -> int:
+        ...
+    def __index__(self) -> int:
+        ...
+    def __init__(self, value: typing.SupportsInt) -> None:
+        ...
+    def __int__(self) -> int:
+        ...
+    def __ne__(self, other: typing.Any) -> bool:
+        ...
+    @typing.overload
+    def __repr__(self) -> str:
+        ...
+    @typing.overload
+    def __repr__(self) -> str:
+        ...
+    def __setstate__(self, state: typing.SupportsInt) -> None:
+        ...
+    def __str__(self) -> str:
+        ...
+    @property
+    def name(self) -> str:
+        ...
+    @property
+    def value(self) -> int:
+        ...
+class StepDerivatives:
+    @property
+    def dYdt(self) -> dict[fourdst._phys.atomic.Species, float]:
+        """
+        The right-hand side (dY/dt) of the ODE system.
+        """
+    @property
+    def energy(self) -> float:
+        """
+        The energy generation rate.
+        """
+def build_nuclear_network(composition: ..., weakInterpolator: ..., maxLayers: gridfire._gridfire.engine.NetworkBuildDepth | typing.SupportsInt = ..., ReactionTypes: NetworkConstructionFlags = ...) -> gridfire._gridfire.reaction.ReactionSet:
+    """
+    Build a nuclear network from a composition using all archived reaction data.
+    """
+def primeNetwork(netIn: gridfire._gridfire.type.NetIn, engine: ..., ignoredReactionTypes: collections.abc.Sequence[...] | None = None) -> PrimingReport:
+    """
+    Prime a network with a short timescale ignition
+    """
+def regularize_jacobian(jacobian: NetworkJacobian, composition: fourdst._phys.composition.Composition) -> NetworkJacobian:
+    """
+    regularize_jacobian
+    """
+ACTIVE: SpeciesStatus  # value = <SpeciesStatus.ACTIVE: 0>
+ADAPTIVE_ENGINE_VIEW: EngineTypes  # value = <EngineTypes.ADAPTIVE_ENGINE_VIEW: 1>
+DEFAULT: NetworkConstructionFlags  # value = <NetworkConstructionFlags.REACLIB: 33>
+DEFINED_ENGINE_VIEW: EngineTypes  # value = <EngineTypes.DEFINED_ENGINE_VIEW: 4>
+EQUILIBRIUM: SpeciesStatus  # value = <SpeciesStatus.EQUILIBRIUM: 1>
+FILE_DEFINED_ENGINE_VIEW: EngineTypes  # value = <EngineTypes.FILE_DEFINED_ENGINE_VIEW: 5>
+FULL_SUCCESS: PrimingReportStatus  # value = <PrimingReportStatus.FULL_SUCCESS: 0>
+FifthOrder: NetworkBuildDepth  # value = <NetworkBuildDepth.FifthOrder: 5>
+FourthOrder: NetworkBuildDepth  # value = <NetworkBuildDepth.FourthOrder: 4>
+Full: NetworkBuildDepth  # value = <NetworkBuildDepth.Full: -1>
+GRAPH_ENGINE: EngineTypes  # value = <EngineTypes.GRAPH_ENGINE: 0>
+INACTIVE_FLOW: SpeciesStatus  # value = <SpeciesStatus.INACTIVE_FLOW: 2>
+MAX_ITERATIONS_REACHED: PrimingReportStatus  # value = <PrimingReportStatus.MAX_ITERATIONS_REACHED: 1>
+MULTISCALE_PARTITIONING_ENGINE_VIEW: EngineTypes  # value = <EngineTypes.MULTISCALE_PARTITIONING_ENGINE_VIEW: 2>
+NONE: NetworkConstructionFlags  # value = <NetworkConstructionFlags.NONE: 0>
+NOT_PRESENT: SpeciesStatus  # value = <SpeciesStatus.NOT_PRESENT: 3>
+NO_SPECIES_TO_PRIME: PrimingReportStatus  # value = <PrimingReportStatus.NO_SPECIES_TO_PRIME: 2>
+PRIMING_ENGINE_VIEW: EngineTypes  # value = <EngineTypes.PRIMING_ENGINE_VIEW: 3>
+REACLIB: NetworkConstructionFlags  # value = <NetworkConstructionFlags.REACLIB: 33>
+REACLIB_STRONG: NetworkConstructionFlags  # value = <NetworkConstructionFlags.REACLIB_STRONG: 1>
+REACLIB_WEAK: NetworkConstructionFlags  # value = <NetworkConstructionFlags.REACLIB_WEAK: 32>
+SecondOrder: NetworkBuildDepth  # value = <NetworkBuildDepth.SecondOrder: 2>
+Shallow: NetworkBuildDepth  # value = <NetworkBuildDepth.Shallow: 1>
+ThirdOrder: NetworkBuildDepth  # value = <NetworkBuildDepth.ThirdOrder: 3>
+WRL_BETA_MINUS: NetworkConstructionFlags  # value = <NetworkConstructionFlags.WRL_BETA_MINUS: 2>
+WRL_BETA_PLUS: NetworkConstructionFlags  # value = <NetworkConstructionFlags.WRL_BETA_PLUS: 4>
+WRL_ELECTRON_CAPTURE: NetworkConstructionFlags  # value = <NetworkConstructionFlags.WRL_ELECTRON_CAPTURE: 8>
+WRL_POSITRON_CAPTURE: NetworkConstructionFlags  # value = <NetworkConstructionFlags.WRL_POSITRON_CAPTURE: 16>
+WRL_WEAK: NetworkConstructionFlags  # value = <NetworkConstructionFlags.WRL_WEAK: 30>
diff --git a/stubs/gridfire/_gridfire/engine/diagnostics.pyi b/stubs/gridfire/_gridfire/engine/diagnostics.pyi
new file mode 100644
index 00000000..5e25cf40
--- /dev/null
+++ b/stubs/gridfire/_gridfire/engine/diagnostics.pyi
@@ -0,0 +1,15 @@
+"""
+A submodule for engine diagnostics
+"""
+from __future__ import annotations
+import collections.abc
+import fourdst._phys.composition
+import gridfire._gridfire.engine
+import typing
+__all__: list[str] = ['inspect_jacobian_stiffness', 'inspect_species_balance', 'report_limiting_species']
+def inspect_jacobian_stiffness(engine: gridfire._gridfire.engine.DynamicEngine, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat, json: bool) -> ... | None:
+    ...
+def inspect_species_balance(engine: gridfire._gridfire.engine.DynamicEngine, species_name: str, comp: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat, json: bool) -> ... | None:
+    ...
+def report_limiting_species(engine: gridfire._gridfire.engine.DynamicEngine, Y_full: collections.abc.Sequence[typing.SupportsFloat], E_full: collections.abc.Sequence[typing.SupportsFloat], relTol: typing.SupportsFloat, absTol: typing.SupportsFloat, top_n: typing.SupportsInt, json: bool) -> ... | None:
+    ...
diff --git a/stubs/gridfire/_gridfire/exceptions.pyi b/stubs/gridfire/_gridfire/exceptions.pyi
new file mode 100644
index 00000000..bfcb1ea1
--- /dev/null
+++ b/stubs/gridfire/_gridfire/exceptions.pyi
@@ -0,0 +1,59 @@
+"""
+GridFire exceptions bindings
+"""
+from __future__ import annotations
+__all__: list[str] = ['BadCollectionError', 'BadRHSEngineError', 'CVODESolverFailureError', 'DebugException', 'EngineError', 'FailedToPartitionEngineError', 'GridFireError', 'HashingError', 'IllConditionedJacobianError', 'InvalidQSESolutionError', 'JacobianError', 'KINSolSolverFailureError', 'MissingBaseReactionError', 'MissingKeyReactionError', 'MissingSeedSpeciesError', 'NetworkResizedError', 'PolicyError', 'ReactionError', 'ReactionParsingError', 'SUNDIALSError', 'SingularJacobianError', 'SolverError', 'StaleJacobianError', 'UnableToSetNetworkReactionsError', 'UninitializedJacobianError', 'UnknownJacobianError', 'UtilityError']
+class BadCollectionError(EngineError):
+    pass
+class BadRHSEngineError(EngineError):
+    pass
+class CVODESolverFailureError(SUNDIALSError):
+    pass
+class DebugException(GridFireError):
+    pass
+class EngineError(GridFireError):
+    pass
+class FailedToPartitionEngineError(EngineError):
+    pass
+class GridFireError(Exception):
+    pass
+class HashingError(UtilityError):
+    pass
+class IllConditionedJacobianError(SolverError):
+    pass
+class InvalidQSESolutionError(EngineError):
+    pass
+class JacobianError(EngineError):
+    pass
+class KINSolSolverFailureError(SUNDIALSError):
+    pass
+class MissingBaseReactionError(PolicyError):
+    pass
+class MissingKeyReactionError(PolicyError):
+    pass
+class MissingSeedSpeciesError(PolicyError):
+    pass
+class NetworkResizedError(EngineError):
+    pass
+class PolicyError(GridFireError):
+    pass
+class ReactionError(GridFireError):
+    pass
+class ReactionParsingError(ReactionError):
+    pass
+class SUNDIALSError(SolverError):
+    pass
+class SingularJacobianError(SolverError):
+    pass
+class SolverError(GridFireError):
+    pass
+class StaleJacobianError(JacobianError):
+    pass
+class UnableToSetNetworkReactionsError(EngineError):
+    pass
+class UninitializedJacobianError(JacobianError):
+    pass
+class UnknownJacobianError(JacobianError):
+    pass
+class UtilityError(GridFireError):
+    pass
diff --git a/stubs/gridfire/_gridfire/io.pyi b/stubs/gridfire/_gridfire/io.pyi
new file mode 100644
index 00000000..6024a26b
--- /dev/null
+++ b/stubs/gridfire/_gridfire/io.pyi
@@ -0,0 +1,14 @@
+"""
+GridFire io bindings
+"""
+from __future__ import annotations
+__all__: list[str] = ['NetworkFileParser', 'ParsedNetworkData', 'SimpleReactionListFileParser']
+class NetworkFileParser:
+    pass
+class ParsedNetworkData:
+    pass
+class SimpleReactionListFileParser(NetworkFileParser):
+    def parse(self, filename: str) -> ParsedNetworkData:
+        """
+        Parse a simple reaction list file and return a ParsedNetworkData object.
+        """
diff --git a/stubs/gridfire/_gridfire/partition.pyi b/stubs/gridfire/_gridfire/partition.pyi
new file mode 100644
index 00000000..50c78f24
--- /dev/null
+++ b/stubs/gridfire/_gridfire/partition.pyi
@@ -0,0 +1,142 @@
+"""
+GridFire partition function bindings
+"""
+from __future__ import annotations
+import collections.abc
+import typing
+__all__: list[str] = ['BasePartitionType', 'CompositePartitionFunction', 'GroundState', 'GroundStatePartitionFunction', 'PartitionFunction', 'RauscherThielemann', 'RauscherThielemannPartitionDataRecord', 'RauscherThielemannPartitionFunction', 'basePartitionTypeToString', 'stringToBasePartitionType']
+class BasePartitionType:
+    """
+    Members:
+    
+      RauscherThielemann
+    
+      GroundState
+    """
+    GroundState: typing.ClassVar[BasePartitionType]  # value = <BasePartitionType.GroundState: 1>
+    RauscherThielemann: typing.ClassVar[BasePartitionType]  # value = <BasePartitionType.RauscherThielemann: 0>
+    __members__: typing.ClassVar[dict[str, BasePartitionType]]  # value = {'RauscherThielemann': <BasePartitionType.RauscherThielemann: 0>, 'GroundState': <BasePartitionType.GroundState: 1>}
+    def __eq__(self, other: typing.Any) -> bool:
+        ...
+    def __getstate__(self) -> int:
+        ...
+    def __hash__(self) -> int:
+        ...
+    def __index__(self) -> int:
+        ...
+    def __init__(self, value: typing.SupportsInt) -> None:
+        ...
+    def __int__(self) -> int:
+        ...
+    def __ne__(self, other: typing.Any) -> bool:
+        ...
+    def __repr__(self) -> str:
+        ...
+    def __setstate__(self, state: typing.SupportsInt) -> None:
+        ...
+    def __str__(self) -> str:
+        ...
+    @property
+    def name(self) -> str:
+        ...
+    @property
+    def value(self) -> int:
+        ...
+class CompositePartitionFunction:
+    @typing.overload
+    def __init__(self, partitionFunctions: collections.abc.Sequence[BasePartitionType]) -> None:
+        """
+        Create a composite partition function from a list of base partition types.
+        """
+    @typing.overload
+    def __init__(self, arg0: CompositePartitionFunction) -> None:
+        """
+        Copy constructor for CompositePartitionFunction.
+        """
+    def evaluate(self, z: typing.SupportsInt, a: typing.SupportsInt, T9: typing.SupportsFloat) -> float:
+        """
+        Evaluate the composite partition function for given Z, A, and T9.
+        """
+    def evaluateDerivative(self, z: typing.SupportsInt, a: typing.SupportsInt, T9: typing.SupportsFloat) -> float:
+        """
+        Evaluate the derivative of the composite partition function for given Z, A, and T9.
+        """
+    def get_type(self) -> str:
+        """
+        Get the type of the partition function (should return 'Composite').
+        """
+    def supports(self, z: typing.SupportsInt, a: typing.SupportsInt) -> bool:
+        """
+        Check if the composite partition function supports given Z and A.
+        """
+class GroundStatePartitionFunction(PartitionFunction):
+    def __init__(self) -> None:
+        ...
+    def evaluate(self, z: typing.SupportsInt, a: typing.SupportsInt, T9: typing.SupportsFloat) -> float:
+        """
+        Evaluate the ground state partition function for given Z, A, and T9.
+        """
+    def evaluateDerivative(self, z: typing.SupportsInt, a: typing.SupportsInt, T9: typing.SupportsFloat) -> float:
+        """
+        Evaluate the derivative of the ground state partition function for given Z, A, and T9.
+        """
+    def get_type(self) -> str:
+        """
+        Get the type of the partition function (should return 'GroundState').
+        """
+    def supports(self, z: typing.SupportsInt, a: typing.SupportsInt) -> bool:
+        """
+        Check if the ground state partition function supports given Z and A.
+        """
+class PartitionFunction:
+    pass
+class RauscherThielemannPartitionDataRecord:
+    @property
+    def a(self) -> int:
+        """
+        Mass number
+        """
+    @property
+    def ground_state_spin(self) -> float:
+        """
+        Ground state spin
+        """
+    @property
+    def normalized_g_values(self) -> float:
+        """
+        Normalized g-values for the first 24 energy levels
+        """
+    @property
+    def z(self) -> int:
+        """
+        Atomic number
+        """
+class RauscherThielemannPartitionFunction(PartitionFunction):
+    def __init__(self) -> None:
+        ...
+    def evaluate(self, z: typing.SupportsInt, a: typing.SupportsInt, T9: typing.SupportsFloat) -> float:
+        """
+        Evaluate the Rauscher-Thielemann partition function for given Z, A, and T9.
+        """
+    def evaluateDerivative(self, z: typing.SupportsInt, a: typing.SupportsInt, T9: typing.SupportsFloat) -> float:
+        """
+        Evaluate the derivative of the Rauscher-Thielemann partition function for given Z, A, and T9.
+        """
+    def get_type(self) -> str:
+        """
+        Get the type of the partition function (should return 'RauscherThielemann').
+        """
+    def supports(self, z: typing.SupportsInt, a: typing.SupportsInt) -> bool:
+        """
+        Check if the Rauscher-Thielemann partition function supports given Z and A.
+        """
+def basePartitionTypeToString(type: BasePartitionType) -> str:
+    """
+    Convert BasePartitionType to string.
+    """
+def stringToBasePartitionType(typeStr: str) -> BasePartitionType:
+    """
+    Convert string to BasePartitionType.
+    """
+GroundState: BasePartitionType  # value = <BasePartitionType.GroundState: 1>
+RauscherThielemann: BasePartitionType  # value = <BasePartitionType.RauscherThielemann: 0>
diff --git a/stubs/gridfire/_gridfire/policy.pyi b/stubs/gridfire/_gridfire/policy.pyi
new file mode 100644
index 00000000..b5bf1af1
--- /dev/null
+++ b/stubs/gridfire/_gridfire/policy.pyi
@@ -0,0 +1,750 @@
+"""
+GridFire network policy bindings
+"""
+from __future__ import annotations
+import collections.abc
+import fourdst._phys.atomic
+import fourdst._phys.composition
+import gridfire._gridfire.engine
+import gridfire._gridfire.reaction
+import typing
+__all__: list[str] = ['CNOChainPolicy', 'CNOIChainPolicy', 'CNOIIChainPolicy', 'CNOIIIChainPolicy', 'CNOIVChainPolicy', 'HotCNOChainPolicy', 'HotCNOIChainPolicy', 'HotCNOIIChainPolicy', 'HotCNOIIIChainPolicy', 'INITIALIZED_UNVERIFIED', 'INITIALIZED_VERIFIED', 'MISSING_KEY_REACTION', 'MISSING_KEY_SPECIES', 'MainSequencePolicy', 'MainSequenceReactionChainPolicy', 'MultiReactionChainPolicy', 'NetworkPolicy', 'NetworkPolicyStatus', 'ProtonProtonChainPolicy', 'ProtonProtonIChainPolicy', 'ProtonProtonIIChainPolicy', 'ProtonProtonIIIChainPolicy', 'ReactionChainPolicy', 'TemperatureDependentChainPolicy', 'TripleAlphaChainPolicy', 'UNINITIALIZED']
+class CNOChainPolicy(MultiReactionChainPolicy):
+    def __eq__(self, other: ReactionChainPolicy) -> bool:
+        """
+        Check equality with another ReactionChainPolicy.
+        """
+    def __hash__(self) -> int:
+        ...
+    def __init__(self) -> None:
+        ...
+    def __ne__(self, other: ReactionChainPolicy) -> bool:
+        """
+        Check inequality with another ReactionChainPolicy.
+        """
+    def __repr__(self) -> str:
+        ...
+    @typing.overload
+    def contains(self, id: str) -> bool:
+        """
+        Check if the reaction chain contains a reaction with the given ID.
+        """
+    @typing.overload
+    def contains(self, reaction: ...) -> bool:
+        """
+        Check if the reaction chain contains the given reaction.
+        """
+    def get_reactions(self) -> gridfire._gridfire.reaction.ReactionSet:
+        """
+        Get the ReactionSet representing this reaction chain.
+        """
+    def hash(self, seed: typing.SupportsInt) -> int:
+        """
+        Compute a hash value for the reaction chain policy.
+        """
+    @typing.overload
+    def name(self) -> str:
+        """
+        Get the name of the reaction chain policy.
+        """
+    @typing.overload
+    def name(self) -> str:
+        """
+        Get the name of the reaction chain policy.
+        """
+class CNOIChainPolicy(TemperatureDependentChainPolicy):
+    def __eq__(self, other: ReactionChainPolicy) -> bool:
+        """
+        Check equality with another ReactionChainPolicy.
+        """
+    def __hash__(self) -> int:
+        ...
+    def __init__(self) -> None:
+        ...
+    def __ne__(self, other: ReactionChainPolicy) -> bool:
+        """
+        Check inequality with another ReactionChainPolicy.
+        """
+    def __repr__(self) -> str:
+        ...
+    @typing.overload
+    def contains(self, id: str) -> bool:
+        """
+        Check if the reaction chain contains a reaction with the given ID.
+        """
+    @typing.overload
+    def contains(self, reaction: ...) -> bool:
+        """
+        Check if the reaction chain contains the given reaction.
+        """
+    def get_reactions(self) -> gridfire._gridfire.reaction.ReactionSet:
+        """
+        Get the ReactionSet representing this reaction chain.
+        """
+    def hash(self, seed: typing.SupportsInt) -> int:
+        """
+        Compute a hash value for the reaction chain policy.
+        """
+    @typing.overload
+    def name(self) -> str:
+        """
+        Get the name of the reaction chain policy.
+        """
+    @typing.overload
+    def name(self) -> str:
+        """
+        Get the name of the reaction chain policy.
+        """
+class CNOIIChainPolicy(TemperatureDependentChainPolicy):
+    def __eq__(self, other: ReactionChainPolicy) -> bool:
+        """
+        Check equality with another ReactionChainPolicy.
+        """
+    def __hash__(self) -> int:
+        ...
+    def __init__(self) -> None:
+        ...
+    def __ne__(self, other: ReactionChainPolicy) -> bool:
+        """
+        Check inequality with another ReactionChainPolicy.
+        """
+    def __repr__(self) -> str:
+        ...
+    @typing.overload
+    def contains(self, id: str) -> bool:
+        """
+        Check if the reaction chain contains a reaction with the given ID.
+        """
+    @typing.overload
+    def contains(self, reaction: ...) -> bool:
+        """
+        Check if the reaction chain contains the given reaction.
+        """
+    def get_reactions(self) -> gridfire._gridfire.reaction.ReactionSet:
+        """
+        Get the ReactionSet representing this reaction chain.
+        """
+    def hash(self, seed: typing.SupportsInt) -> int:
+        """
+        Compute a hash value for the reaction chain policy.
+        """
+    @typing.overload
+    def name(self) -> str:
+        """
+        Get the name of the reaction chain policy.
+        """
+    @typing.overload
+    def name(self) -> str:
+        """
+        Get the name of the reaction chain policy.
+        """
+class CNOIIIChainPolicy(TemperatureDependentChainPolicy):
+    def __eq__(self, other: ReactionChainPolicy) -> bool:
+        """
+        Check equality with another ReactionChainPolicy.
+        """
+    def __hash__(self) -> int:
+        ...
+    def __init__(self) -> None:
+        ...
+    def __ne__(self, other: ReactionChainPolicy) -> bool:
+        """
+        Check inequality with another ReactionChainPolicy.
+        """
+    def __repr__(self) -> str:
+        ...
+    @typing.overload
+    def contains(self, id: str) -> bool:
+        """
+        Check if the reaction chain contains a reaction with the given ID.
+        """
+    @typing.overload
+    def contains(self, reaction: ...) -> bool:
+        """
+        Check if the reaction chain contains the given reaction.
+        """
+    def get_reactions(self) -> gridfire._gridfire.reaction.ReactionSet:
+        """
+        Get the ReactionSet representing this reaction chain.
+        """
+    def hash(self, seed: typing.SupportsInt) -> int:
+        """
+        Compute a hash value for the reaction chain policy.
+        """
+    @typing.overload
+    def name(self) -> str:
+        """
+        Get the name of the reaction chain policy.
+        """
+    @typing.overload
+    def name(self) -> str:
+        """
+        Get the name of the reaction chain policy.
+        """
+class CNOIVChainPolicy(TemperatureDependentChainPolicy):
+    def __eq__(self, other: ReactionChainPolicy) -> bool:
+        """
+        Check equality with another ReactionChainPolicy.
+        """
+    def __hash__(self) -> int:
+        ...
+    def __init__(self) -> None:
+        ...
+    def __ne__(self, other: ReactionChainPolicy) -> bool:
+        """
+        Check inequality with another ReactionChainPolicy.
+        """
+    def __repr__(self) -> str:
+        ...
+    @typing.overload
+    def contains(self, id: str) -> bool:
+        """
+        Check if the reaction chain contains a reaction with the given ID.
+        """
+    @typing.overload
+    def contains(self, reaction: ...) -> bool:
+        """
+        Check if the reaction chain contains the given reaction.
+        """
+    def get_reactions(self) -> gridfire._gridfire.reaction.ReactionSet:
+        """
+        Get the ReactionSet representing this reaction chain.
+        """
+    def hash(self, seed: typing.SupportsInt) -> int:
+        """
+        Compute a hash value for the reaction chain policy.
+        """
+    @typing.overload
+    def name(self) -> str:
+        """
+        Get the name of the reaction chain policy.
+        """
+    @typing.overload
+    def name(self) -> str:
+        """
+        Get the name of the reaction chain policy.
+        """
+class HotCNOChainPolicy(MultiReactionChainPolicy):
+    def __eq__(self, other: ReactionChainPolicy) -> bool:
+        """
+        Check equality with another ReactionChainPolicy.
+        """
+    def __hash__(self) -> int:
+        ...
+    def __init__(self) -> None:
+        ...
+    def __ne__(self, other: ReactionChainPolicy) -> bool:
+        """
+        Check inequality with another ReactionChainPolicy.
+        """
+    def __repr__(self) -> str:
+        ...
+    @typing.overload
+    def contains(self, id: str) -> bool:
+        """
+        Check if the reaction chain contains a reaction with the given ID.
+        """
+    @typing.overload
+    def contains(self, reaction: ...) -> bool:
+        """
+        Check if the reaction chain contains the given reaction.
+        """
+    def get_reactions(self) -> gridfire._gridfire.reaction.ReactionSet:
+        """
+        Get the ReactionSet representing this reaction chain.
+        """
+    def hash(self, seed: typing.SupportsInt) -> int:
+        """
+        Compute a hash value for the reaction chain policy.
+        """
+    @typing.overload
+    def name(self) -> str:
+        """
+        Get the name of the reaction chain policy.
+        """
+    @typing.overload
+    def name(self) -> str:
+        """
+        Get the name of the reaction chain policy.
+        """
+class HotCNOIChainPolicy(TemperatureDependentChainPolicy):
+    def __eq__(self, other: ReactionChainPolicy) -> bool:
+        """
+        Check equality with another ReactionChainPolicy.
+        """
+    def __hash__(self) -> int:
+        ...
+    def __init__(self) -> None:
+        ...
+    def __ne__(self, other: ReactionChainPolicy) -> bool:
+        """
+        Check inequality with another ReactionChainPolicy.
+        """
+    def __repr__(self) -> str:
+        ...
+    @typing.overload
+    def contains(self, id: str) -> bool:
+        """
+        Check if the reaction chain contains a reaction with the given ID.
+        """
+    @typing.overload
+    def contains(self, reaction: ...) -> bool:
+        """
+        Check if the reaction chain contains the given reaction.
+        """
+    def get_reactions(self) -> gridfire._gridfire.reaction.ReactionSet:
+        """
+        Get the ReactionSet representing this reaction chain.
+        """
+    def hash(self, seed: typing.SupportsInt) -> int:
+        """
+        Compute a hash value for the reaction chain policy.
+        """
+    @typing.overload
+    def name(self) -> str:
+        """
+        Get the name of the reaction chain policy.
+        """
+    @typing.overload
+    def name(self) -> str:
+        """
+        Get the name of the reaction chain policy.
+        """
+class HotCNOIIChainPolicy(TemperatureDependentChainPolicy):
+    def __eq__(self, other: ReactionChainPolicy) -> bool:
+        """
+        Check equality with another ReactionChainPolicy.
+        """
+    def __hash__(self) -> int:
+        ...
+    def __init__(self) -> None:
+        ...
+    def __ne__(self, other: ReactionChainPolicy) -> bool:
+        """
+        Check inequality with another ReactionChainPolicy.
+        """
+    def __repr__(self) -> str:
+        ...
+    @typing.overload
+    def contains(self, id: str) -> bool:
+        """
+        Check if the reaction chain contains a reaction with the given ID.
+        """
+    @typing.overload
+    def contains(self, reaction: ...) -> bool:
+        """
+        Check if the reaction chain contains the given reaction.
+        """
+    def get_reactions(self) -> gridfire._gridfire.reaction.ReactionSet:
+        """
+        Get the ReactionSet representing this reaction chain.
+        """
+    def hash(self, seed: typing.SupportsInt) -> int:
+        """
+        Compute a hash value for the reaction chain policy.
+        """
+    @typing.overload
+    def name(self) -> str:
+        """
+        Get the name of the reaction chain policy.
+        """
+    @typing.overload
+    def name(self) -> str:
+        """
+        Get the name of the reaction chain policy.
+        """
+class HotCNOIIIChainPolicy(TemperatureDependentChainPolicy):
+    def __eq__(self, other: ReactionChainPolicy) -> bool:
+        """
+        Check equality with another ReactionChainPolicy.
+        """
+    def __hash__(self) -> int:
+        ...
+    def __init__(self) -> None:
+        ...
+    def __ne__(self, other: ReactionChainPolicy) -> bool:
+        """
+        Check inequality with another ReactionChainPolicy.
+        """
+    def __repr__(self) -> str:
+        ...
+    @typing.overload
+    def contains(self, id: str) -> bool:
+        """
+        Check if the reaction chain contains a reaction with the given ID.
+        """
+    @typing.overload
+    def contains(self, reaction: ...) -> bool:
+        """
+        Check if the reaction chain contains the given reaction.
+        """
+    def get_reactions(self) -> gridfire._gridfire.reaction.ReactionSet:
+        """
+        Get the ReactionSet representing this reaction chain.
+        """
+    def hash(self, seed: typing.SupportsInt) -> int:
+        """
+        Compute a hash value for the reaction chain policy.
+        """
+    @typing.overload
+    def name(self) -> str:
+        """
+        Get the name of the reaction chain policy.
+        """
+    @typing.overload
+    def name(self) -> str:
+        """
+        Get the name of the reaction chain policy.
+        """
+class MainSequencePolicy(NetworkPolicy):
+    @typing.overload
+    def __init__(self, composition: fourdst._phys.composition.Composition) -> None:
+        """
+        Construct MainSequencePolicy from an existing composition.
+        """
+    @typing.overload
+    def __init__(self, seed_species: collections.abc.Sequence[fourdst._phys.atomic.Species], mass_fractions: collections.abc.Sequence[typing.SupportsFloat]) -> None:
+        """
+        Construct MainSequencePolicy from seed species and mass fractions.
+        """
+    def construct(self) -> gridfire._gridfire.engine.DynamicEngine:
+        """
+        Construct the network according to the policy.
+        """
+    def get_engine_types_stack(self) -> list[gridfire._gridfire.engine.EngineTypes]:
+        """
+        Get the types of engines in the stack constructed by the network policy.
+        """
+    def get_seed_reactions(self) -> gridfire._gridfire.reaction.ReactionSet:
+        """
+        Get the set of seed reactions required by the network policy.
+        """
+    def get_seed_species(self) -> set[fourdst._phys.atomic.Species]:
+        """
+        Get the set of seed species required by the network policy.
+        """
+    def get_status(self) -> NetworkPolicyStatus:
+        """
+        Get the current status of the network policy.
+        """
+    def name(self) -> str:
+        """
+        Get the name of the network policy.
+        """
+class MainSequenceReactionChainPolicy(MultiReactionChainPolicy):
+    def __eq__(self, other: ReactionChainPolicy) -> bool:
+        """
+        Check equality with another ReactionChainPolicy.
+        """
+    def __hash__(self) -> int:
+        ...
+    def __init__(self) -> None:
+        ...
+    def __ne__(self, other: ReactionChainPolicy) -> bool:
+        """
+        Check inequality with another ReactionChainPolicy.
+        """
+    def __repr__(self) -> str:
+        ...
+    @typing.overload
+    def contains(self, id: str) -> bool:
+        """
+        Check if the reaction chain contains a reaction with the given ID.
+        """
+    @typing.overload
+    def contains(self, reaction: ...) -> bool:
+        """
+        Check if the reaction chain contains the given reaction.
+        """
+    def get_reactions(self) -> gridfire._gridfire.reaction.ReactionSet:
+        """
+        Get the ReactionSet representing this reaction chain.
+        """
+    def hash(self, seed: typing.SupportsInt) -> int:
+        """
+        Compute a hash value for the reaction chain policy.
+        """
+    @typing.overload
+    def name(self) -> str:
+        """
+        Get the name of the reaction chain policy.
+        """
+    @typing.overload
+    def name(self) -> str:
+        """
+        Get the name of the reaction chain policy.
+        """
+class MultiReactionChainPolicy(ReactionChainPolicy):
+    pass
+class NetworkPolicy:
+    pass
+class NetworkPolicyStatus:
+    """
+    Members:
+    
+      UNINITIALIZED
+    
+      INITIALIZED_UNVERIFIED
+    
+      MISSING_KEY_REACTION
+    
+      MISSING_KEY_SPECIES
+    
+      INITIALIZED_VERIFIED
+    """
+    INITIALIZED_UNVERIFIED: typing.ClassVar[NetworkPolicyStatus]  # value = <NetworkPolicyStatus.INITIALIZED_UNVERIFIED: 1>
+    INITIALIZED_VERIFIED: typing.ClassVar[NetworkPolicyStatus]  # value = <NetworkPolicyStatus.INITIALIZED_VERIFIED: 4>
+    MISSING_KEY_REACTION: typing.ClassVar[NetworkPolicyStatus]  # value = <NetworkPolicyStatus.MISSING_KEY_REACTION: 2>
+    MISSING_KEY_SPECIES: typing.ClassVar[NetworkPolicyStatus]  # value = <NetworkPolicyStatus.MISSING_KEY_SPECIES: 3>
+    UNINITIALIZED: typing.ClassVar[NetworkPolicyStatus]  # value = <NetworkPolicyStatus.UNINITIALIZED: 0>
+    __members__: typing.ClassVar[dict[str, NetworkPolicyStatus]]  # value = {'UNINITIALIZED': <NetworkPolicyStatus.UNINITIALIZED: 0>, 'INITIALIZED_UNVERIFIED': <NetworkPolicyStatus.INITIALIZED_UNVERIFIED: 1>, 'MISSING_KEY_REACTION': <NetworkPolicyStatus.MISSING_KEY_REACTION: 2>, 'MISSING_KEY_SPECIES': <NetworkPolicyStatus.MISSING_KEY_SPECIES: 3>, 'INITIALIZED_VERIFIED': <NetworkPolicyStatus.INITIALIZED_VERIFIED: 4>}
+    def __eq__(self, other: typing.Any) -> bool:
+        ...
+    def __getstate__(self) -> int:
+        ...
+    def __hash__(self) -> int:
+        ...
+    def __index__(self) -> int:
+        ...
+    def __init__(self, value: typing.SupportsInt) -> None:
+        ...
+    def __int__(self) -> int:
+        ...
+    def __ne__(self, other: typing.Any) -> bool:
+        ...
+    def __repr__(self) -> str:
+        ...
+    def __setstate__(self, state: typing.SupportsInt) -> None:
+        ...
+    def __str__(self) -> str:
+        ...
+    @property
+    def name(self) -> str:
+        ...
+    @property
+    def value(self) -> int:
+        ...
+class ProtonProtonChainPolicy(MultiReactionChainPolicy):
+    def __eq__(self, other: ReactionChainPolicy) -> bool:
+        """
+        Check equality with another ReactionChainPolicy.
+        """
+    def __hash__(self) -> int:
+        ...
+    def __init__(self) -> None:
+        ...
+    def __ne__(self, other: ReactionChainPolicy) -> bool:
+        """
+        Check inequality with another ReactionChainPolicy.
+        """
+    def __repr__(self) -> str:
+        ...
+    @typing.overload
+    def contains(self, id: str) -> bool:
+        """
+        Check if the reaction chain contains a reaction with the given ID.
+        """
+    @typing.overload
+    def contains(self, reaction: ...) -> bool:
+        """
+        Check if the reaction chain contains the given reaction.
+        """
+    def get_reactions(self) -> gridfire._gridfire.reaction.ReactionSet:
+        """
+        Get the ReactionSet representing this reaction chain.
+        """
+    def hash(self, seed: typing.SupportsInt) -> int:
+        """
+        Compute a hash value for the reaction chain policy.
+        """
+    @typing.overload
+    def name(self) -> str:
+        """
+        Get the name of the reaction chain policy.
+        """
+    @typing.overload
+    def name(self) -> str:
+        """
+        Get the name of the reaction chain policy.
+        """
+class ProtonProtonIChainPolicy(TemperatureDependentChainPolicy):
+    def __eq__(self, other: ReactionChainPolicy) -> bool:
+        """
+        Check equality with another ReactionChainPolicy.
+        """
+    def __hash__(self) -> int:
+        ...
+    def __init__(self) -> None:
+        ...
+    def __ne__(self, other: ReactionChainPolicy) -> bool:
+        """
+        Check inequality with another ReactionChainPolicy.
+        """
+    def __repr__(self) -> str:
+        ...
+    @typing.overload
+    def contains(self, id: str) -> bool:
+        """
+        Check if the reaction chain contains a reaction with the given ID.
+        """
+    @typing.overload
+    def contains(self, reaction: ...) -> bool:
+        """
+        Check if the reaction chain contains the given reaction.
+        """
+    def get_reactions(self) -> gridfire._gridfire.reaction.ReactionSet:
+        """
+        Get the ReactionSet representing this reaction chain.
+        """
+    def hash(self, seed: typing.SupportsInt) -> int:
+        """
+        Compute a hash value for the reaction chain policy.
+        """
+    @typing.overload
+    def name(self) -> str:
+        """
+        Get the name of the reaction chain policy.
+        """
+    @typing.overload
+    def name(self) -> str:
+        """
+        Get the name of the reaction chain policy.
+        """
+class ProtonProtonIIChainPolicy(TemperatureDependentChainPolicy):
+    def __eq__(self, other: ReactionChainPolicy) -> bool:
+        """
+        Check equality with another ReactionChainPolicy.
+        """
+    def __hash__(self) -> int:
+        ...
+    def __init__(self) -> None:
+        ...
+    def __ne__(self, other: ReactionChainPolicy) -> bool:
+        """
+        Check inequality with another ReactionChainPolicy.
+        """
+    def __repr__(self) -> str:
+        ...
+    @typing.overload
+    def contains(self, id: str) -> bool:
+        """
+        Check if the reaction chain contains a reaction with the given ID.
+        """
+    @typing.overload
+    def contains(self, reaction: ...) -> bool:
+        """
+        Check if the reaction chain contains the given reaction.
+        """
+    def get_reactions(self) -> gridfire._gridfire.reaction.ReactionSet:
+        """
+        Get the ReactionSet representing this reaction chain.
+        """
+    def hash(self, seed: typing.SupportsInt) -> int:
+        """
+        Compute a hash value for the reaction chain policy.
+        """
+    @typing.overload
+    def name(self) -> str:
+        """
+        Get the name of the reaction chain policy.
+        """
+    @typing.overload
+    def name(self) -> str:
+        """
+        Get the name of the reaction chain policy.
+        """
+class ProtonProtonIIIChainPolicy(TemperatureDependentChainPolicy):
+    def __eq__(self, other: ReactionChainPolicy) -> bool:
+        """
+        Check equality with another ReactionChainPolicy.
+        """
+    def __hash__(self) -> int:
+        ...
+    def __init__(self) -> None:
+        ...
+    def __ne__(self, other: ReactionChainPolicy) -> bool:
+        """
+        Check inequality with another ReactionChainPolicy.
+        """
+    def __repr__(self) -> str:
+        ...
+    @typing.overload
+    def contains(self, id: str) -> bool:
+        """
+        Check if the reaction chain contains a reaction with the given ID.
+        """
+    @typing.overload
+    def contains(self, reaction: ...) -> bool:
+        """
+        Check if the reaction chain contains the given reaction.
+        """
+    def get_reactions(self) -> gridfire._gridfire.reaction.ReactionSet:
+        """
+        Get the ReactionSet representing this reaction chain.
+        """
+    def hash(self, seed: typing.SupportsInt) -> int:
+        """
+        Compute a hash value for the reaction chain policy.
+        """
+    @typing.overload
+    def name(self) -> str:
+        """
+        Get the name of the reaction chain policy.
+        """
+    @typing.overload
+    def name(self) -> str:
+        """
+        Get the name of the reaction chain policy.
+        """
+class ReactionChainPolicy:
+    pass
+class TemperatureDependentChainPolicy(ReactionChainPolicy):
+    pass
+class TripleAlphaChainPolicy(TemperatureDependentChainPolicy):
+    def __eq__(self, other: ReactionChainPolicy) -> bool:
+        """
+        Check equality with another ReactionChainPolicy.
+        """
+    def __hash__(self) -> int:
+        ...
+    def __init__(self) -> None:
+        ...
+    def __ne__(self, other: ReactionChainPolicy) -> bool:
+        """
+        Check inequality with another ReactionChainPolicy.
+        """
+    def __repr__(self) -> str:
+        ...
+    @typing.overload
+    def contains(self, id: str) -> bool:
+        """
+        Check if the reaction chain contains a reaction with the given ID.
+        """
+    @typing.overload
+    def contains(self, reaction: ...) -> bool:
+        """
+        Check if the reaction chain contains the given reaction.
+        """
+    def get_reactions(self) -> gridfire._gridfire.reaction.ReactionSet:
+        """
+        Get the ReactionSet representing this reaction chain.
+        """
+    def hash(self, seed: typing.SupportsInt) -> int:
+        """
+        Compute a hash value for the reaction chain policy.
+        """
+    @typing.overload
+    def name(self) -> str:
+        """
+        Get the name of the reaction chain policy.
+        """
+    @typing.overload
+    def name(self) -> str:
+        """
+        Get the name of the reaction chain policy.
+        """
+INITIALIZED_UNVERIFIED: NetworkPolicyStatus  # value = <NetworkPolicyStatus.INITIALIZED_UNVERIFIED: 1>
+INITIALIZED_VERIFIED: NetworkPolicyStatus  # value = <NetworkPolicyStatus.INITIALIZED_VERIFIED: 4>
+MISSING_KEY_REACTION: NetworkPolicyStatus  # value = <NetworkPolicyStatus.MISSING_KEY_REACTION: 2>
+MISSING_KEY_SPECIES: NetworkPolicyStatus  # value = <NetworkPolicyStatus.MISSING_KEY_SPECIES: 3>
+UNINITIALIZED: NetworkPolicyStatus  # value = <NetworkPolicyStatus.UNINITIALIZED: 0>
diff --git a/stubs/gridfire/_gridfire/reaction.pyi b/stubs/gridfire/_gridfire/reaction.pyi
new file mode 100644
index 00000000..d2659168
--- /dev/null
+++ b/stubs/gridfire/_gridfire/reaction.pyi
@@ -0,0 +1,249 @@
+"""
+GridFire reaction bindings
+"""
+from __future__ import annotations
+import collections.abc
+import fourdst._phys.atomic
+import fourdst._phys.composition
+import typing
+__all__: list[str] = ['LogicalReaclibReaction', 'RateCoefficientSet', 'ReaclibReaction', 'ReactionSet', 'get_all_reactions', 'packReactionSet']
+class LogicalReaclibReaction(ReaclibReaction):
+    @typing.overload
+    def __init__(self, reactions: collections.abc.Sequence[ReaclibReaction]) -> None:
+        """
+        Construct a LogicalReaclibReaction from a vector of ReaclibReaction objects.
+        """
+    @typing.overload
+    def __init__(self, reactions: collections.abc.Sequence[ReaclibReaction], is_reverse: bool) -> None:
+        """
+        Construct a LogicalReaclibReaction from a vector of ReaclibReaction objects.
+        """
+    def __len__(self) -> int:
+        """
+        Overload len() to return the number of source rates.
+        """
+    def add_reaction(self, reaction: ReaclibReaction) -> None:
+        """
+        Add another Reaction source to this logical reaction.
+        """
+    def calculate_forward_rate_log_derivative(self, T9: typing.SupportsFloat, rho: typing.SupportsFloat, Ye: typing.SupportsFloat, mue: typing.SupportsFloat, Composition: fourdst._phys.composition.Composition) -> float:
+        """
+        Calculate the forward rate log derivative at a given temperature T9 (in units of 10^9 K).
+        """
+    def calculate_rate(self, T9: typing.SupportsFloat, rho: typing.SupportsFloat, Ye: typing.SupportsFloat, mue: typing.SupportsFloat, Y: collections.abc.Sequence[typing.SupportsFloat], index_to_species_map: collections.abc.Mapping[typing.SupportsInt, fourdst._phys.atomic.Species]) -> float:
+        """
+        Calculate the reaction rate at a given temperature T9 (in units of 10^9 K). Note that for a reaclib reaction only T9 is actually used, all other parameters are there for interface compatibility.
+        """
+    def size(self) -> int:
+        """
+        Get the number of source rates contributing to this logical reaction.
+        """
+    def sources(self) -> list[str]:
+        """
+        Get the list of source labels for the aggregated rates.
+        """
+class RateCoefficientSet:
+    def __init__(self, a0: typing.SupportsFloat, a1: typing.SupportsFloat, a2: typing.SupportsFloat, a3: typing.SupportsFloat, a4: typing.SupportsFloat, a5: typing.SupportsFloat, a6: typing.SupportsFloat) -> None:
+        """
+        Construct a RateCoefficientSet with the given parameters.
+        """
+class ReaclibReaction:
+    __hash__: typing.ClassVar[None] = None
+    def __eq__(self, arg0: ReaclibReaction) -> bool:
+        """
+        Equality operator for reactions based on their IDs.
+        """
+    def __init__(self, id: str, peName: str, chapter: typing.SupportsInt, reactants: collections.abc.Sequence[fourdst._phys.atomic.Species], products: collections.abc.Sequence[fourdst._phys.atomic.Species], qValue: typing.SupportsFloat, label: str, sets: RateCoefficientSet, reverse: bool = False) -> None:
+        """
+        Construct a Reaction with the given parameters.
+        """
+    def __neq__(self, arg0: ReaclibReaction) -> bool:
+        """
+        Inequality operator for reactions based on their IDs.
+        """
+    def __repr__(self) -> str:
+        ...
+    def all_species(self) -> set[fourdst._phys.atomic.Species]:
+        """
+        Get all species involved in the reaction (both reactants and products) as a set.
+        """
+    def calculate_rate(self, T9: typing.SupportsFloat, rho: typing.SupportsFloat, Y: collections.abc.Sequence[typing.SupportsFloat]) -> float:
+        """
+        Calculate the reaction rate at a given temperature T9 (in units of 10^9 K).
+        """
+    def chapter(self) -> int:
+        """
+        Get the REACLIB chapter number defining the reaction structure.
+        """
+    def contains(self, species: fourdst._phys.atomic.Species) -> bool:
+        """
+        Check if the reaction contains a specific species.
+        """
+    def contains_product(self, arg0: fourdst._phys.atomic.Species) -> bool:
+        """
+        Check if the reaction contains a specific product species.
+        """
+    def contains_reactant(self, arg0: fourdst._phys.atomic.Species) -> bool:
+        """
+        Check if the reaction contains a specific reactant species.
+        """
+    def excess_energy(self) -> float:
+        """
+        Calculate the excess energy from the mass difference of reactants and products.
+        """
+    def hash(self, seed: typing.SupportsInt = 0) -> int:
+        """
+        Compute a hash for the reaction based on its ID.
+        """
+    def id(self) -> str:
+        """
+        Get the unique identifier of the reaction.
+        """
+    def is_reverse(self) -> bool:
+        """
+        Check if this is a reverse reaction rate.
+        """
+    def num_species(self) -> int:
+        """
+        Count the number of species in the reaction.
+        """
+    def peName(self) -> str:
+        """
+        Get the reaction name in (projectile, ejectile) notation (e.g., 'p(p,g)d').
+        """
+    def product_species(self) -> set[fourdst._phys.atomic.Species]:
+        """
+        Get the product species of the reaction as a set.
+        """
+    def products(self) -> list[fourdst._phys.atomic.Species]:
+        """
+        Get a list of product species in the reaction.
+        """
+    def qValue(self) -> float:
+        """
+        Get the Q-value of the reaction in MeV.
+        """
+    def rateCoefficients(self) -> RateCoefficientSet:
+        """
+        get the set of rate coefficients.
+        """
+    def reactant_species(self) -> set[fourdst._phys.atomic.Species]:
+        """
+        Get the reactant species of the reaction as a set.
+        """
+    def reactants(self) -> list[fourdst._phys.atomic.Species]:
+        """
+        Get a list of reactant species in the reaction.
+        """
+    def sourceLabel(self) -> str:
+        """
+        Get the source label for the rate data (e.g., 'wc12w', 'st08').
+        """
+    @typing.overload
+    def stoichiometry(self, species: fourdst._phys.atomic.Species) -> int:
+        """
+        Get the stoichiometry of the reaction as a map from species to their coefficients.
+        """
+    @typing.overload
+    def stoichiometry(self) -> dict[fourdst._phys.atomic.Species, int]:
+        """
+        Get the stoichiometry of the reaction as a map from species to their coefficients.
+        """
+class ReactionSet:
+    __hash__: typing.ClassVar[None] = None
+    @staticmethod
+    def from_clones(reactions: collections.abc.Sequence[...]) -> ReactionSet:
+        """
+        Create a ReactionSet that takes ownership of the reactions by cloning the input reactions.
+        """
+    def __eq__(self, LogicalReactionSet: ReactionSet) -> bool:
+        """
+        Equality operator for LogicalReactionSets based on their contents.
+        """
+    def __getitem__(self, index: typing.SupportsInt) -> ...:
+        """
+        Get a LogicalReaclibReaction by index.
+        """
+    def __getitem___(self, id: str) -> ...:
+        """
+        Get a LogicalReaclibReaction by its ID.
+        """
+    @typing.overload
+    def __init__(self, reactions: collections.abc.Sequence[...]) -> None:
+        """
+        Construct a LogicalReactionSet from a vector of LogicalReaclibReaction objects.
+        """
+    @typing.overload
+    def __init__(self) -> None:
+        """
+        Default constructor for an empty LogicalReactionSet.
+        """
+    @typing.overload
+    def __init__(self, other: ReactionSet) -> None:
+        """
+        Copy constructor for LogicalReactionSet.
+        """
+    def __len__(self) -> int:
+        """
+        Overload len() to return the number of LogicalReactions.
+        """
+    def __ne__(self, LogicalReactionSet: ReactionSet) -> bool:
+        """
+        Inequality operator for LogicalReactionSets based on their contents.
+        """
+    def __repr__(self) -> str:
+        ...
+    def add_reaction(self, reaction: ...) -> None:
+        """
+        Add a LogicalReaclibReaction to the set.
+        """
+    def clear(self) -> None:
+        """
+        Remove all LogicalReactions from the set.
+        """
+    @typing.overload
+    def contains(self, id: str) -> bool:
+        """
+        Check if the set contains a specific LogicalReaclibReaction.
+        """
+    @typing.overload
+    def contains(self, reaction: ...) -> bool:
+        """
+        Check if the set contains a specific Reaction.
+        """
+    def contains_product(self, species: fourdst._phys.atomic.Species) -> bool:
+        """
+        Check if any reaction in the set has the species as a product.
+        """
+    def contains_reactant(self, species: fourdst._phys.atomic.Species) -> bool:
+        """
+        Check if any reaction in the set has the species as a reactant.
+        """
+    def contains_species(self, species: fourdst._phys.atomic.Species) -> bool:
+        """
+        Check if any reaction in the set involves the given species.
+        """
+    def getReactionSetSpecies(self) -> set[fourdst._phys.atomic.Species]:
+        """
+        Get all species involved in the reactions of the set as a set of Species objects.
+        """
+    def hash(self, seed: typing.SupportsInt = 0) -> int:
+        """
+        Compute a hash for the LogicalReactionSet based on its contents.
+        """
+    def remove_reaction(self, reaction: ...) -> None:
+        """
+        Remove a LogicalReaclibReaction from the set.
+        """
+    def size(self) -> int:
+        """
+        Get the number of LogicalReactions in the set.
+        """
+def get_all_reactions() -> ReactionSet:
+    """
+    Get all reactions from the REACLIB database.
+    """
+def packReactionSet(reactionSet: ReactionSet) -> ReactionSet:
+    """
+    Convert a ReactionSet to a LogicalReactionSet by aggregating reactions with the same peName.
+    """
diff --git a/stubs/gridfire/_gridfire/screening.pyi b/stubs/gridfire/_gridfire/screening.pyi
new file mode 100644
index 00000000..e2b9c181
--- /dev/null
+++ b/stubs/gridfire/_gridfire/screening.pyi
@@ -0,0 +1,68 @@
+"""
+GridFire plasma screening bindings
+"""
+from __future__ import annotations
+import collections.abc
+import fourdst._phys.atomic
+import gridfire._gridfire.reaction
+import typing
+__all__: list[str] = ['BARE', 'BareScreeningModel', 'ScreeningModel', 'ScreeningType', 'WEAK', 'WeakScreeningModel', 'selectScreeningModel']
+class BareScreeningModel:
+    def __init__(self) -> None:
+        ...
+    def calculateScreeningFactors(self, reactions: gridfire._gridfire.reaction.ReactionSet, species: collections.abc.Sequence[fourdst._phys.atomic.Species], Y: collections.abc.Sequence[typing.SupportsFloat], T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> list[float]:
+        """
+        Calculate the bare plasma screening factors. This always returns 1.0 (bare)
+        """
+class ScreeningModel:
+    pass
+class ScreeningType:
+    """
+    Members:
+    
+      BARE
+    
+      WEAK
+    """
+    BARE: typing.ClassVar[ScreeningType]  # value = <ScreeningType.BARE: 0>
+    WEAK: typing.ClassVar[ScreeningType]  # value = <ScreeningType.WEAK: 1>
+    __members__: typing.ClassVar[dict[str, ScreeningType]]  # value = {'BARE': <ScreeningType.BARE: 0>, 'WEAK': <ScreeningType.WEAK: 1>}
+    def __eq__(self, other: typing.Any) -> bool:
+        ...
+    def __getstate__(self) -> int:
+        ...
+    def __hash__(self) -> int:
+        ...
+    def __index__(self) -> int:
+        ...
+    def __init__(self, value: typing.SupportsInt) -> None:
+        ...
+    def __int__(self) -> int:
+        ...
+    def __ne__(self, other: typing.Any) -> bool:
+        ...
+    def __repr__(self) -> str:
+        ...
+    def __setstate__(self, state: typing.SupportsInt) -> None:
+        ...
+    def __str__(self) -> str:
+        ...
+    @property
+    def name(self) -> str:
+        ...
+    @property
+    def value(self) -> int:
+        ...
+class WeakScreeningModel:
+    def __init__(self) -> None:
+        ...
+    def calculateScreeningFactors(self, reactions: gridfire._gridfire.reaction.ReactionSet, species: collections.abc.Sequence[fourdst._phys.atomic.Species], Y: collections.abc.Sequence[typing.SupportsFloat], T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> list[float]:
+        """
+        Calculate the weak plasma screening factors using the Salpeter (1954) model.
+        """
+def selectScreeningModel(type: ScreeningType) -> ScreeningModel:
+    """
+    Select a screening model based on the specified type. Returns a pointer to the selected model.
+    """
+BARE: ScreeningType  # value = <ScreeningType.BARE: 0>
+WEAK: ScreeningType  # value = <ScreeningType.WEAK: 1>
diff --git a/stubs/gridfire/_gridfire/solver.pyi b/stubs/gridfire/_gridfire/solver.pyi
new file mode 100644
index 00000000..ee9938ef
--- /dev/null
+++ b/stubs/gridfire/_gridfire/solver.pyi
@@ -0,0 +1,92 @@
+"""
+GridFire numerical solver bindings
+"""
+from __future__ import annotations
+import collections.abc
+import fourdst._phys.atomic
+import gridfire._gridfire.engine
+import gridfire._gridfire.type
+import typing
+__all__: list[str] = ['CVODESolverStrategy', 'CVODETimestepContext', 'DynamicNetworkSolverStrategy', 'SolverContextBase']
+class CVODESolverStrategy(DynamicNetworkSolverStrategy):
+    def __init__(self, engine: gridfire._gridfire.engine.DynamicEngine) -> None:
+        """
+        Initialize the CVODESolverStrategy object.
+        """
+    def evaluate(self, netIn: gridfire._gridfire.type.NetIn, display_trigger: bool = False) -> gridfire._gridfire.type.NetOut:
+        """
+        evaluate the dynamic engine using the dynamic engine class
+        """
+    def get_absTol(self) -> float:
+        """
+        Get the absolute tolerance for the CVODE solver.
+        """
+    def get_relTol(self) -> float:
+        """
+        Get the relative tolerance for the CVODE solver.
+        """
+    def get_stdout_logging_enabled(self) -> bool:
+        """
+        Check if solver logging to standard output is enabled.
+        """
+    def set_absTol(self, absTol: typing.SupportsFloat) -> None:
+        """
+        Set the absolute tolerance for the CVODE solver.
+        """
+    def set_callback(self, cb: collections.abc.Callable[[CVODETimestepContext], None]) -> None:
+        """
+        Set a callback function which will run at the end of every successful timestep
+        """
+    def set_relTol(self, relTol: typing.SupportsFloat) -> None:
+        """
+        Set the relative tolerance for the CVODE solver.
+        """
+    def set_stdout_logging_enabled(self, logging_enabled: bool) -> None:
+        """
+        Enable logging to standard output.
+        """
+class CVODETimestepContext(SolverContextBase):
+    @property
+    def T9(self) -> float:
+        ...
+    @property
+    def currentConvergenceFailures(self) -> int:
+        ...
+    @property
+    def currentNonlinearIterations(self) -> int:
+        ...
+    @property
+    def dt(self) -> float:
+        ...
+    @property
+    def engine(self) -> gridfire._gridfire.engine.DynamicEngine:
+        ...
+    @property
+    def last_step_time(self) -> float:
+        ...
+    @property
+    def networkSpecies(self) -> list[fourdst._phys.atomic.Species]:
+        ...
+    @property
+    def num_steps(self) -> int:
+        ...
+    @property
+    def rho(self) -> float:
+        ...
+    @property
+    def state(self) -> list[float]:
+        ...
+    @property
+    def t(self) -> float:
+        ...
+class DynamicNetworkSolverStrategy:
+    def describe_callback_context(self) -> list[tuple[str, str]]:
+        """
+        Get a structure representing what data is in the callback context in a human readable format
+        """
+    def evaluate(self, netIn: gridfire._gridfire.type.NetIn) -> gridfire._gridfire.type.NetOut:
+        """
+        evaluate the dynamic engine using the dynamic engine class
+        """
+class SolverContextBase:
+    pass
diff --git a/stubs/gridfire/_gridfire/type.pyi b/stubs/gridfire/_gridfire/type.pyi
new file mode 100644
index 00000000..45332f62
--- /dev/null
+++ b/stubs/gridfire/_gridfire/type.pyi
@@ -0,0 +1,61 @@
+"""
+GridFire type bindings
+"""
+from __future__ import annotations
+import fourdst._phys.composition
+import typing
+__all__: list[str] = ['NetIn', 'NetOut']
+class NetIn:
+    composition: fourdst._phys.composition.Composition
+    def __init__(self) -> None:
+        ...
+    def __repr__(self) -> str:
+        ...
+    @property
+    def density(self) -> float:
+        ...
+    @density.setter
+    def density(self, arg0: typing.SupportsFloat) -> None:
+        ...
+    @property
+    def dt0(self) -> float:
+        ...
+    @dt0.setter
+    def dt0(self, arg0: typing.SupportsFloat) -> None:
+        ...
+    @property
+    def energy(self) -> float:
+        ...
+    @energy.setter
+    def energy(self, arg0: typing.SupportsFloat) -> None:
+        ...
+    @property
+    def tMax(self) -> float:
+        ...
+    @tMax.setter
+    def tMax(self, arg0: typing.SupportsFloat) -> None:
+        ...
+    @property
+    def temperature(self) -> float:
+        ...
+    @temperature.setter
+    def temperature(self, arg0: typing.SupportsFloat) -> None:
+        ...
+class NetOut:
+    def __repr__(self) -> str:
+        ...
+    @property
+    def composition(self) -> fourdst._phys.composition.Composition:
+        ...
+    @property
+    def dEps_dRho(self) -> float:
+        ...
+    @property
+    def dEps_dT(self) -> float:
+        ...
+    @property
+    def energy(self) -> float:
+        ...
+    @property
+    def num_steps(self) -> int:
+        ...
diff --git a/stubs/gridfire/_gridfire/utils/__init__.pyi b/stubs/gridfire/_gridfire/utils/__init__.pyi
new file mode 100644
index 00000000..d8738d61
--- /dev/null
+++ b/stubs/gridfire/_gridfire/utils/__init__.pyi
@@ -0,0 +1,17 @@
+"""
+GridFire utility method bindings
+"""
+from __future__ import annotations
+import fourdst._phys.composition
+import gridfire._gridfire.engine
+import typing
+from . import hashing
+__all__: list[str] = ['formatNuclearTimescaleLogString', 'hash_atomic', 'hash_reaction', 'hashing']
+def formatNuclearTimescaleLogString(engine: gridfire._gridfire.engine.DynamicEngine, Y: fourdst._phys.composition.Composition, T9: typing.SupportsFloat, rho: typing.SupportsFloat) -> str:
+    """
+    Format a string for logging nuclear timescales based on temperature, density, and energy generation rate.
+    """
+def hash_atomic(a: typing.SupportsInt, z: typing.SupportsInt) -> int:
+    ...
+def hash_reaction(reaction: ...) -> int:
+    ...
diff --git a/stubs/gridfire/_gridfire/utils/hashing/__init__.pyi b/stubs/gridfire/_gridfire/utils/hashing/__init__.pyi
new file mode 100644
index 00000000..b85a400a
--- /dev/null
+++ b/stubs/gridfire/_gridfire/utils/hashing/__init__.pyi
@@ -0,0 +1,6 @@
+"""
+module for gridfire hashing functions
+"""
+from __future__ import annotations
+from . import reaction
+__all__: list[str] = ['reaction']
diff --git a/stubs/gridfire/_gridfire/utils/hashing/reaction.pyi b/stubs/gridfire/_gridfire/utils/hashing/reaction.pyi
new file mode 100644
index 00000000..131f482d
--- /dev/null
+++ b/stubs/gridfire/_gridfire/utils/hashing/reaction.pyi
@@ -0,0 +1,12 @@
+"""
+utility module for hashing gridfire reaction functions
+"""
+from __future__ import annotations
+import typing
+__all__: list[str] = ['mix_species', 'multiset_combine', 'splitmix64']
+def mix_species(a: typing.SupportsInt, z: typing.SupportsInt) -> int:
+    ...
+def multiset_combine(acc: typing.SupportsInt, x: typing.SupportsInt) -> int:
+    ...
+def splitmix64(x: typing.SupportsInt) -> int:
+    ...
diff --git a/subprojects/fourdst.wrap b/subprojects/fourdst.wrap
index 4cd10a0a..7048ec12 100644
--- a/subprojects/fourdst.wrap
+++ b/subprojects/fourdst.wrap
@@ -1,4 +1,4 @@
 [wrap-git]
 url = https://github.com/4D-STAR/fourdst
-revision = v0.9.5
+revision = v0.9.6
 depth = 1
diff --git a/tests/graphnet_sandbox/main.cpp b/tests/graphnet_sandbox/main.cpp
index 3171f9cd..aaa04783 100644
--- a/tests/graphnet_sandbox/main.cpp
+++ b/tests/graphnet_sandbox/main.cpp
@@ -270,6 +270,10 @@ int main(int argc, char** argv) {
     CLI11_PARSE(app, argc, argv);
 
     const NetIn netIn = init(temp, rho, tMax);
+    std::println("Starting Integration with T = {} K, ρ = {} g/cm³, tMax = {} s", temp, rho, tMax);
+    std::println("Composition is: {}", utils::iterable_to_delimited_string(netIn.composition, ", ", [&netIn](std::pair<const fourdst::atomic::Species&, double> arg) {
+        return std::format("{:5}: {:10.4E}", arg.first.name(), arg.second);
+    }));
 
     policy::MainSequencePolicy stellarPolicy(netIn.composition);
     stellarPolicy.construct();
diff --git a/validation/pynucastro/GridFireEquiv/GridFireEvolve.py b/validation/pynucastro/GridFireEquiv/GridFireEvolve.py
new file mode 100644
index 00000000..1a5db34d
--- /dev/null
+++ b/validation/pynucastro/GridFireEquiv/GridFireEvolve.py
@@ -0,0 +1,95 @@
+from fourdst.composition import Composition
+from gridfire.type import NetIn
+from gridfire.policy import MainSequencePolicy
+from gridfire.solver import CVODESolverStrategy
+from enum import Enum
+from typing import Dict, Union, SupportsFloat
+import json
+import dicttoxml
+
+def init_composition() -> Composition:
+    Y = [7.0262E-01, 9.7479E-06, 6.8955E-02, 2.5000E-04, 7.8554E-05, 6.0144E-04, 8.1031E-05, 2.1513E-05]
+    S = ["H-1", "He-3", "He-4", "C-12", "N-14", "O-16", "Ne-20", "Mg-24"]
+    return Composition(S, Y)
+
+def init_netIn(temp: float, rho: float, time: float, comp: Composition) -> NetIn:
+    netIn = NetIn()
+    netIn.temperature = temp
+    netIn.density = rho
+    netIn.tMax = time
+    netIn.dt0 = 1e-12
+    netIn.composition = comp
+    return netIn
+
+class StepData(Enum):
+    TIME = 0
+    DT = 1
+    COMP = 2
+    CONTRIB = 3
+
+
+class StepLogger:
+    def __init__(self):
+        self.num_steps: int = 0
+        self.step_data: Dict[int, Dict[StepData, Union[SupportsFloat, Dict[str, SupportsFloat]]]] = {}
+
+    def log_step(self, context):
+        engine = context.engine
+        self.step_data[self.num_steps] = {}
+        self.step_data[self.num_steps][StepData.TIME] = context.t
+        self.step_data[self.num_steps][StepData.DT] = context.dt
+        comp_data: Dict[str, SupportsFloat] = {}
+        for species in engine.getNetworkSpecies():
+            sid = engine.getSpeciesIndex(species)
+            comp_data[species.name()] = context.state[sid]
+        self.step_data[self.num_steps][StepData.COMP] = comp_data
+        self.num_steps += 1
+
+    def to_json (self, filename: str):
+        serializable_data = {
+            stepNum: {
+                StepData.TIME.name: step[StepData.TIME],
+                StepData.DT.name: step[StepData.DT],
+                StepData.COMP.name: step[StepData.COMP],
+            }
+            for stepNum, step in self.step_data.items()
+        }
+        with open(filename, 'w') as f:
+            json.dump(serializable_data, f, indent=4)
+
+    def to_xml(self, filename: str):
+        serializable_data = {
+            stepNum: {
+                StepData.TIME.name: step[StepData.TIME],
+                StepData.DT.name: step[StepData.DT],
+                StepData.COMP.name: step[StepData.COMP],
+            }
+            for stepNum, step in self.step_data.items()
+        }
+        xml_data = dicttoxml.dicttoxml(serializable_data, custom_root='StepLog', attr_type=False)
+        with open(filename, 'wb') as f:
+            f.write(xml_data)
+
+def main(temp: float, rho: float, time: float):
+    comp = init_composition()
+    netIn = init_netIn(temp, rho, time, comp)
+
+    policy = MainSequencePolicy(comp)
+    engine = policy.construct()
+
+    solver = CVODESolverStrategy(engine)
+
+    step_logger = StepLogger()
+    solver.set_callback(lambda context: step_logger.log_step(context))
+
+    solver.evaluate(netIn, False)
+    step_logger.to_xml("log_data.xml")
+
+if __name__ == "__main__":
+    import argparse
+    parser = argparse.ArgumentParser(description="Simple python example of GridFire usage")
+    parser.add_argument("-t", "--temp", type=float, help="Temperature in K", default=1.5e7)
+    parser.add_argument("-r", "--rho", type=float, help="Density in g/cm^3", default=1.5e2)
+    parser.add_argument("--tMax", type=float, help="Time in s", default=3.1536 * 1e17)
+    args = parser.parse_args()
+    main(args.temp, args.rho, args.tMax)