feat(io): added ability to export the rate calculations of any engine to a python file

GridFire can now (to a very limited degree) approximate the abilities of pynucastro in so far as it can not just solve an engine but it can also export it to a python file. Currently only the rates are exported so no solver can yet be used (since things like the rest of the RHS and the jacobian are missing) but this is a step in that direction

src/include/gridfire/io/generative/python.h

@@ -0,0 +1,21 @@
+#pragma once
+
+#include <string>
+#include <vector>
+
+#include "gridfire/reaction/reaction.h"
+#include "gridfire/engine/engine_abstract.h"
+
+namespace gridfire::io::gen {
+    struct PyFunctionDef {
+        std::string func_name;
+        std::string func_code;
+        std::vector<std::string> module_req;
+    };
+
+    PyFunctionDef exportReactionToPy(const reaction::Reaction& reaction);
+
+    std::string exportEngineToPy(const DynamicEngine& engine);
+
+    void exportEngineToPy(const DynamicEngine& engine, const std::string& fileName);
+}

src/include/gridfire/reaction/reaction.h

@@ -344,6 +344,9 @@ namespace gridfire::reaction {
             os << "Reaction(ID: " << r.id() << ")";
             return os;
         }
+
+        virtual std::optional<std::vector<RateCoefficientSet>> getRateCoefficients() const = 0;
+
     };
     class ReaclibReaction : public Reaction {
     public:
@@ -443,6 +446,8 @@ namespace gridfire::reaction {
          */
         [[nodiscard]] const RateCoefficientSet& rateCoefficients() const { return m_rateCoefficients; }
 
+        [[nodiscard]] std::optional<std::vector<RateCoefficientSet>> getRateCoefficients() const override;
+
         /**
          * @brief Checks if the reaction involves a given species as a reactant or product.
          * @param species The species to check for.
@@ -695,6 +700,8 @@ namespace gridfire::reaction {
             CppAD::AD<double> mue, const std::vector<CppAD::AD<double>>& Y, const std::unordered_map<size_t,fourdst::atomic::Species>& index_to_species_map
         ) const override;
 
+        [[nodiscard]] std::optional<std::vector<RateCoefficientSet>> getRateCoefficients() const override;
+
         /** @name Iterators
          *  Provides iterators to loop over the rate coefficient sets.
          */

src/include/gridfire/reaction/weak/weak.h

@@ -361,6 +361,8 @@ namespace gridfire::rates::weak {
          */
         [[nodiscard]] const WeakRateInterpolator& getWeakRateInterpolator() const;
 
+        [[nodiscard]] std::optional<std::vector<reaction::RateCoefficientSet>> getRateCoefficients() const override { return std::nullopt; }
+
     private:
         /**
          * @brief Internal unified implementation for scalar/AD rate evaluation.

src/lib/io/generative/python.cpp

@@ -0,0 +1,162 @@
+#include "gridfire/io/generative/python.h"
+#include "fourdst/composition/species.h"
+
+#include <string>
+#include <vector>
+#include <format>
+#include <algorithm>
+#include <ranges>
+#include <fstream>
+
+#include <optional>
+
+#include "gridfire/engine/engine_abstract.h"
+
+namespace {
+    template <typename T>
+    std::string join(std::vector<T> arr, std::string delim) {
+        if (arr.empty()) return {};
+
+        size_t total = delim.size() * (arr.size() - 1);
+        for (const auto& e : arr) total += std::string(e).size();
+
+        std::string out;
+        out.reserve(total);
+
+        for (size_t i = 0; i < arr.size(); i++) {
+            out += std::string(arr[i]);
+            if (i < arr.size() - 1) {
+                out += delim;
+            }
+        }
+        return out;
+    }
+
+    std::string format_reaction_function_name(const gridfire::reaction::Reaction& reaction) {
+        const std::vector<fourdst::atomic::Species>& reactants = reaction.reactants();
+        const std::vector<fourdst::atomic::Species>& products = reaction.products();
+
+        auto format_species_names = [](const std::vector<fourdst::atomic::Species>& arr) -> std::vector<std::string> {
+            std::vector<std::string> out;
+            for (const auto& r : arr) {
+                auto name = std::string(r.name());
+                std::erase_if(name, [](const char c) {
+                    return c == '-';
+                });
+                out.push_back(name);
+            }
+            return out;
+        };
+
+        const std::vector<std::string> reactant_names = format_species_names(reactants);
+        const std::vector<std::string> product_names = format_species_names(products);
+
+        std::string funcName = std::format(
+            "{}_to_{}",
+            join<std::string>(reactant_names, "_"),
+            join<std::string>(product_names, "_")
+        );
+
+        return funcName;
+    }
+
+    std::string tfFunc = R"(def tf(T9: float) -> Dict[str, float]:
+    """
+    Calculate T9^(1/3) & T9^(5/3) & ln(T9)
+    """
+    T913 = T9 ** (1.0 / 3.0)
+    T953 = T9 ** (5.0 / 3.0)
+    logT9 = np.log(T9)
+    return {
+        "T913": T913,
+        "T953": T953,
+        "logT9": logT9
+    }
+)";
+}
+
+namespace gridfire::io::gen {
+    PyFunctionDef exportReactionToPy(const reaction::Reaction& reaction) {
+        std::string funcName = format_reaction_function_name(reaction);
+        std::optional<std::vector<reaction::RateCoefficientSet>> rateSets = reaction.getRateCoefficients();
+
+        std::string signature = std::format(
+        "def {}(T9: float) -> float:",
+            funcName
+        );
+        if (!rateSets || rateSets -> empty()) {
+            const std::string funcCode = std::format(
+                "{}\n#No Reaclib Rates Defined\n    return 0.0",
+                signature,
+                funcName
+            );
+            return {funcName, funcCode, {}};
+        }
+
+        std::string commentHeader = std::format(
+            "    # Calculates the reaction rate for the reaction: {}",
+            reaction.id()
+        );
+
+        std::string callTf = R"(    tfv = tf(T9)
+    T913 = tfv["T913"]
+    T953 = tfv["T953"]
+    logT9 = tfv["logT9"]
+)";
+        std::string initRate = "    rate = 0.0\n";
+
+        std::vector<std::string> rateLines;
+        for (const auto&[a0, a1, a2, a3, a4, a5, a6] : *rateSets) {
+            std::string rateLine = std::format(
+                "    rate += np.exp(({}) + ({} / T9) + ({} / T913) + ({} * T913) + ({} * T9) + ({} * T953) + ({} * logT9))",
+                a0,
+                a1,
+                a2,
+                a3,
+                a4,
+                a5,
+                a6
+            );
+            rateLines.push_back(rateLine);
+        }
+
+        const std::string rate_body = join<std::string>(rateLines, "\n");
+
+        const std::string returnCall = "    return rate";
+
+        const std::vector<std::string> functionLines = {
+            signature,
+            commentHeader,
+            callTf,
+            initRate,
+            rate_body,
+            returnCall
+        };
+
+        const std::string funcCode = join<std::string>(functionLines, "\n");
+
+        return {funcName, funcCode, {"np"}};
+
+    }
+
+    std::string exportEngineToPy(const gridfire::DynamicEngine& engine) {
+        auto reactions = engine.getNetworkReactions();
+        std::vector<std::string> functions;
+        functions.push_back(R"(import numpy as np
+from typing import Dict, List, Tuple, Callable)");
+        functions.push_back(tfFunc);
+        for (const auto& reaction : reactions) {
+            PyFunctionDef funcDef = exportReactionToPy(*reaction);
+            functions.push_back(funcDef.func_code);
+        }
+        return join<std::string>(functions, "\n\n");
+    }
+
+    void exportEngineToPy(const DynamicEngine &engine, const std::string &fileName) {
+        const std::string funcCode = exportEngineToPy(engine);
+        std::ofstream outFile(fileName);
+        outFile << funcCode;
+        outFile.close();
+    }
+
+}

src/lib/reaction/reaction.cpp

@@ -95,6 +95,12 @@ namespace gridfire::reaction {
         return d_log_k_fwd_dT9; // Return the derivative of the log rate with respect to T9
     }
 
+    std::optional<std::vector<RateCoefficientSet>> ReaclibReaction::getRateCoefficients() const {
+        std::vector<RateCoefficientSet> rateCoefficients;
+        rateCoefficients.push_back(m_rateCoefficients);
+        return rateCoefficients;
+    }
+
     bool ReaclibReaction::contains(const Species &species) const {
         return contains_reactant(species) || contains_product(species);
     }
@@ -321,6 +327,14 @@ namespace gridfire::reaction {
         return calculate_rate<CppAD::AD<double>>(T9);
     }
 
+    std::optional<std::vector<RateCoefficientSet>> LogicalReaclibReaction::getRateCoefficients() const {
+        std::vector<RateCoefficientSet> rates;
+        for (const auto& rate : m_rates) {
+            rates.push_back(rate);
+        }
+        return rates;
+    }
+
     ReactionSet::ReactionSet(
         std::vector<std::unique_ptr<Reaction>>&& reactions
     ) :

src/meson.build

@@ -15,6 +15,7 @@ gridfire_sources = files(
     'lib/reaction/weak/weak.cpp',
     'lib/reaction/weak/weak_interpolator.cpp',
     'lib/io/network_file.cpp',
+    'lib/io/generative/python.cpp',
     'lib/solver/strategies/CVODE_solver_strategy.cpp',
     'lib/solver/strategies/triggers/engine_partitioning_trigger.cpp',
     'lib/screening/screening_types.cpp',