feat(partition-functions): added framework and some concrete partition functions

GroundState partition function, Rauscher&Thielemann partition function, and composite partition function added

src/network/include/gridfire/partition/composite/partition_composite.h

@@ -0,0 +1,43 @@
+#pragma once
+
+#include "gridfire/partition/partition_abstract.h"
+
+#include "fourdst/logging/logging.h"
+
+#include <string>
+#include <unordered_map>
+#include <vector>
+
+#include <memory>
+
+
+namespace gridfire::partition {
+    enum BasePartitionType {
+        RauscherThielemann, ///< Rauscher-Thielemann partition function
+        GroundState,        ///< Ground state partition function
+    };
+
+    inline std::unordered_map<BasePartitionType, std::string> basePartitionTypeToString = {
+        {RauscherThielemann, "RauscherThielemann"},
+        {GroundState, "GroundState"}
+    };
+
+    inline std::unordered_map<std::string, BasePartitionType> stringToBasePartitionType = {
+        {"RauscherThielemann", RauscherThielemann},
+        {"GroundState", GroundState}
+    };
+
+    class CompositePartitionFunction final : public PartitionFunction {
+    public:
+        explicit CompositePartitionFunction(const std::vector<BasePartitionType>& partitionFunctions);
+        double evaluate(int z, int a, double T9) const override;
+        double evaluateDerivative(int z, int a, double T9) const override;
+        bool supports(int z, int a) const override;
+        std::string type() const override;
+    private:
+        quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log");
+        std::vector<std::unique_ptr<PartitionFunction>> m_partitionFunctions; ///< Set of partition functions to use in the composite partition function.
+    private:
+        std::unique_ptr<PartitionFunction> selectPartitionFunction(const BasePartitionType type) const;
+    };
+}

src/network/include/gridfire/partition/partition_abstract.h

@@ -0,0 +1,14 @@
+#pragma once
+
+#include <string>
+
+namespace gridfire::partition {
+    class PartitionFunction {
+    public:
+        virtual ~PartitionFunction() = default;
+        virtual double evaluate(int z, int a, double T9) const = 0;
+        virtual double evaluateDerivative(int z, int a, double T9) const = 0;
+        virtual bool supports(int z, int a) const = 0;
+        virtual std::string type() const = 0;
+    };
+}

src/network/include/gridfire/partition/partition_ground.h

@@ -0,0 +1,38 @@
+#pragma once
+
+#include "gridfire/partition/partition_abstract.h"
+
+#include "fourdst/logging/logging.h"
+
+#include <unordered_map>
+
+#include "quill/Logger.h"
+
+namespace gridfire::partition {
+    class GroundStatePartitionFunction final : public PartitionFunction {
+    public:
+        GroundStatePartitionFunction();
+        double evaluate(
+            const int z,
+            const int a,
+            const double T9
+        ) const override;
+        double evaluateDerivative(
+            const int z,
+            const int a,
+            const double T9
+        ) const override;
+        bool supports(
+            const int z,
+            const int a
+        ) const override;
+        std::string type() const override { return "GroundState"; }
+    private:
+        quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log");
+        std::unordered_map<int, double> m_ground_state_spin;
+        static constexpr int make_key(
+            const int z,
+        const int a);
+    };
+
+}

src/network/include/gridfire/partition/partition_rauscher_thielemann.h

@@ -0,0 +1,56 @@
+#pragma once
+
+#include "gridfire/partition/partition_abstract.h"
+
+#include "fourdst/logging/logging.h"
+
+#include <string>
+#include <vector>
+#include <unordered_map>
+#include <array>
+#include <utility>
+
+namespace gridfire::partition {
+    class RauscherThielemannPartitionFunction final : public PartitionFunction {
+    public:
+        RauscherThielemannPartitionFunction();
+        double evaluate(int z, int a, double T9) const override;
+        double evaluateDerivative(int z, int a, double T9) const override;
+        bool supports(int z, int a) const override;
+        std::string type() const override { return "RauscherThielemann"; }
+    private:
+        enum Bounds {
+            FRONT,
+            BACK,
+            MIDDLE
+        };
+    private:
+        struct IsotopeData {
+            double ground_state_spin;
+            std::array<double, 24> normalized_g_values;
+        };
+        struct InterpolationPoints {
+            double T9_high;
+            double G_norm_high;
+            double T9_low;
+            double G_norm_low;
+        };
+        struct IdentifiedIsotope {
+            Bounds bound;
+            const IsotopeData& data;
+            size_t upperIndex;
+            size_t lowerIndex;
+        };
+    private:
+        quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log");
+        std::unordered_map<int, IsotopeData> m_partitionData;
+    private:
+        static InterpolationPoints get_interpolation_points(
+            const size_t upper_index,
+            const size_t lower_index,
+            const std::array<double, 24>& normalized_g_values
+        );
+        IdentifiedIsotope find(int z, int a, double T9) const;
+        static constexpr int make_key(int z, int a);
+    };
+}

src/network/include/gridfire/partition/rauscher_thielemann_partition_data_record.h

@@ -0,0 +1,15 @@
+#pragma once
+
+#include <cstdint>
+
+namespace gridfire::partition::record {
+#pragma pack(push, 1)
+    struct RauscherThielemannPartitionDataRecord {
+        uint32_t z; ///< Atomic number
+        uint32_t a; ///< Mass number
+        double ground_state_spin; ///< Ground state spin
+        double partition_function; ///< Partition function value
+        double normalized_g_values[24]; ///< Normalized g-values for the first 24 energy levels
+    };
+#pragma pack(pop)
+}

src/network/lib/partition/composite/partition_composite.cpp

@@ -0,0 +1,101 @@
+#include "gridfire/partition/composite/partition_composite.h"
+
+#include <vector>
+#include <set>
+
+#include "gridfire/partition/partition_ground.h"
+#include "gridfire/partition/partition_rauscher_thielemann.h"
+#include "quill/LogMacros.h"
+
+namespace gridfire::partition {
+    CompositePartitionFunction::CompositePartitionFunction(
+        const std::vector<BasePartitionType>& partitionFunctions
+    ) {
+        for (const auto& type : partitionFunctions) {
+            LOG_TRACE_L2(m_logger, "Adding partition function of type: {}", basePartitionTypeToString[type]);
+            m_partitionFunctions.push_back(selectPartitionFunction(type));
+        }
+    }
+
+    double CompositePartitionFunction::evaluate(int z, int a, double T9) const {
+        LOG_TRACE_L1(m_logger, "Evaluating partition function for Z={} A={} T9={}", z, a, T9);
+        for (const auto& partitionFunction : m_partitionFunctions) {
+            if (partitionFunction->supports(z, a)) {
+                LOG_TRACE_L2(m_logger, "Partition function of type {} supports Z={} A={}", partitionFunction->type(), z, a);
+                return partitionFunction->evaluate(z, a, T9);
+            } else {
+                LOG_TRACE_L2(m_logger, "Partition function of type {} does not support Z={} A={}", partitionFunction->type(), z, a);
+            }
+        }
+        LOG_ERROR(
+            m_logger,
+            "No partition function supports Z={} A={} T9={}. Tried: {}",
+            z,
+            a,
+            T9,
+            type()
+        );
+        throw std::runtime_error("No partition function supports the given Z, A, and T9 values.");
+    }
+
+    double CompositePartitionFunction::evaluateDerivative(int z, int a, double T9) const {
+        for (const auto& partitionFunction : m_partitionFunctions) {
+            if (partitionFunction->supports(z, a)) {
+                LOG_TRACE_L2(m_logger, "Evaluating derivative of partition function for Z={} A={} T9={}", z, a, T9);
+                return partitionFunction->evaluateDerivative(z, a, T9);
+            }
+        }
+        LOG_ERROR(
+            m_logger,
+            "No partition function supports Z={} A={} T9={}. Tried: {}",
+            z,
+            a,
+            T9,
+            type()
+        );
+        throw std::runtime_error("No partition function supports the given Z, A, and T9 values.");
+    }
+
+    bool CompositePartitionFunction::supports(int z, int a) const {
+        for (const auto& partitionFunction : m_partitionFunctions) {
+            if (partitionFunction->supports(z, a)) {
+                LOG_TRACE_L2(m_logger, "Partition function supports Z={} A={}", z, a);
+                return true;
+            }
+        }
+        return false;
+    }
+
+    std::string CompositePartitionFunction::type() const {
+        std::stringstream ss;
+        ss << "CompositePartitionFunction(";
+        int count = 0;
+        for (const auto& partitionFunction : m_partitionFunctions) {
+            ss << partitionFunction->type();
+            if (count < m_partitionFunctions.size() - 1) {
+                ss << ", ";
+            }
+            count++;
+        }
+        ss << ")";
+        std::string types = ss.str();
+        return types;
+    }
+
+    std::unique_ptr<PartitionFunction> CompositePartitionFunction::selectPartitionFunction(
+        const BasePartitionType type
+    ) const {
+        switch (type) {
+            case RauscherThielemann: {
+                return std::make_unique<RauscherThielemannPartitionFunction>();
+            }
+            case GroundState: {
+                return std::make_unique<GroundStatePartitionFunction>();
+            }
+            default: {
+                LOG_ERROR(m_logger, "Unknown partition function type");
+                throw std::runtime_error("Unknown partition function type");
+            }
+        }
+    }
+}

src/network/lib/partition/partition_ground.cpp

@@ -0,0 +1,50 @@
+#include "gridfire/partition/partition_ground.h"
+
+#include <ranges>
+
+#include "fourdst/logging/logging.h"
+#include "fourdst/composition/atomicSpecies.h"
+#include "fourdst/composition/species.h"
+
+#include "quill/LogMacros.h"
+
+namespace gridfire::partition {
+    GroundStatePartitionFunction::GroundStatePartitionFunction() {
+        for (const auto &isotope: fourdst::atomic::species | std::views::values) {
+            m_ground_state_spin[make_key(isotope.z(), isotope.a())] = isotope.spin();
+        }
+    }
+
+    double GroundStatePartitionFunction::evaluate(
+        const int z,
+        const int a,
+        const double T9
+    ) const {
+        LOG_TRACE_L2(m_logger, "Evaluating ground state partition function for Z={} A={} T9={}", z, a, T9);
+        const int key = make_key(z, a);
+        const double spin = m_ground_state_spin.at(key);
+        return (2.0 * spin) + 1.0;
+    }
+
+    double GroundStatePartitionFunction::evaluateDerivative(
+        const int z,
+        const int a,
+        const double T9
+    ) const {
+        LOG_TRACE_L2(m_logger, "Evaluating derivative of ground state partition function for Z={} A={} T9={}", z, a, T9);
+        return 0.0;
+    }
+
+    bool GroundStatePartitionFunction::supports(
+        const int z,
+        const int a
+    ) const {
+        return m_ground_state_spin.contains(make_key(z, a));
+    }
+    constexpr int GroundStatePartitionFunction::make_key(
+        const int z,
+        const int a
+    ) {
+        return z * 1000 + a; // Simple key generation for Z and A
+    }
+}

src/network/lib/partition/partition_rauscher_thielemann.cpp

@@ -0,0 +1,141 @@
+#include "gridfire/partition/partition_rauscher_thielemann.h"
+#include "gridfire/partition/rauscher_thielemann_partition_data.h"
+#include "gridfire/partition/rauscher_thielemann_partition_data_record.h"
+
+#include "fourdst/logging/logging.h"
+#include "quill/LogMacros.h"
+
+#include <stdexcept>
+#include <algorithm>
+#include <vector>
+#include <array>
+
+namespace gridfire::partition {
+    static constexpr std::array<double, 24> RT_TEMPERATURE_GRID_T9 = {
+        0.01, 0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.5,
+        2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0
+    };
+
+    RauscherThielemannPartitionFunction::RauscherThielemannPartitionFunction() {
+        constexpr size_t numRecords = rauscher_thielemann_partition_data_len / sizeof(record::RauscherThielemannPartitionDataRecord);
+        m_partitionData.reserve(numRecords);
+        const auto* records = reinterpret_cast<const record::RauscherThielemannPartitionDataRecord*>(rauscher_thielemann_partition_data);
+        for (size_t i = 0; i < numRecords; ++i) {
+            const auto& record = records[i];
+            IsotopeData data;
+            data.ground_state_spin = record.ground_state_spin;
+            std::ranges::copy(record.normalized_g_values, data.normalized_g_values.begin());
+
+            m_partitionData[make_key(record.z, record.a)] = std::move(data);
+        }
+    }
+
+    double RauscherThielemannPartitionFunction::evaluate(
+        const int z,
+        const int a,
+        const double T9
+    ) const {
+        LOG_TRACE_L2(m_logger, "Evaluating Rauscher-Thielemann partition function for Z={} A={} T9={}", z, a, T9);
+
+        const auto [bound, data, upperIndex, lowerIndex] = find(z, a, T9);
+
+        switch (bound) {
+            case FRONT: {
+                LOG_TRACE_L2(m_logger, "Using FRONT bound for Z={} A={} T9={}", z, a, T9);
+                return data.normalized_g_values.front() * (2.0 * data.ground_state_spin + 1.0);
+            }
+            case BACK: {
+                LOG_TRACE_L2(m_logger, "Using BACK bound for Z={} A={} T9={}", z, a, T9);
+                return data.normalized_g_values.back() * (2.0 * data.ground_state_spin + 1.0);
+            }
+            case MIDDLE: {
+                LOG_TRACE_L2(m_logger, "Using MIDDLE bound for Z={} A={} T9={}", z, a, T9);
+            }
+        }
+
+        const auto [T9_high, G_norm_high, T9_low, G_norm_low] = get_interpolation_points(
+            upperIndex,
+            lowerIndex,
+            data.normalized_g_values
+        );
+
+        const double frac = (T9 - T9_low) / (T9_high - T9_low);
+        const double interpolated_g_norm = G_norm_low + frac * (G_norm_high - G_norm_low);
+
+        return interpolated_g_norm * (2.0 * data.ground_state_spin + 1.0);
+    }
+
+    double RauscherThielemannPartitionFunction::evaluateDerivative(
+        const int z,
+        const int a,
+        const double T9
+    ) const {
+        LOG_TRACE_L2(m_logger, "Evaluating derivative of Rauscher-Thielemann partition function for Z={} A={} T9={}", z, a, T9);
+        const auto [bound, data, upperIndex, lowerIndex] = find(z, a, T9);
+        if (bound == FRONT || bound == BACK) {
+            LOG_TRACE_L2(m_logger, "Derivative is zero for Z={} A={} T9={} (bound: {})", z, a, T9, bound == FRONT ? "FRONT" : "BACK");
+            return 0.0; // Derivative is zero at the boundaries
+        }
+        const auto [T9_high, G_norm_high, T9_low, G_norm_low] = get_interpolation_points(
+            upperIndex,
+            lowerIndex,
+            data.normalized_g_values
+        );
+        const double slope_g_norm = (G_norm_high - G_norm_low) / (T9_high - T9_low);
+        return slope_g_norm * (2.0 * data.ground_state_spin + 1.0);
+    }
+
+    bool RauscherThielemannPartitionFunction::supports(
+        const int z,
+        const int a
+    ) const {
+        return m_partitionData.contains(make_key(z, a));
+    }
+
+    RauscherThielemannPartitionFunction::InterpolationPoints RauscherThielemannPartitionFunction::get_interpolation_points(
+        const size_t upper_index,
+        const size_t lower_index,
+        const std::array<double, 24>& normalized_g_values
+    ) {
+        const double T_high = RT_TEMPERATURE_GRID_T9[upper_index];
+        const double G_norm_high = normalized_g_values[upper_index];
+        const double T_low = RT_TEMPERATURE_GRID_T9[lower_index];
+        const double G_norm_low = normalized_g_values[lower_index];
+        return {T_high, G_norm_high, T_low, G_norm_low};
+    }
+
+    RauscherThielemannPartitionFunction::IdentifiedIsotope RauscherThielemannPartitionFunction::find(
+        const int z,
+        const int a,
+        const double T9
+    ) const {
+        const auto key = make_key(z, a);
+        const auto it = m_partitionData.find(key);
+
+        if (it == m_partitionData.end()) {
+            LOG_ERROR(m_logger, "Rauscher-Thielemann partition function data for Z={} A={} not found.", z, a);
+            throw std::out_of_range("Partition function data not found for Z=" + std::to_string(z) + " A=" + std::to_string(a));
+        }
+
+        const IsotopeData& data = it->second;
+        const auto upper_it = std::ranges::lower_bound(RT_TEMPERATURE_GRID_T9, T9);
+        Bounds bound;
+        if (upper_it == RT_TEMPERATURE_GRID_T9.begin()) {
+            bound = FRONT; // T9 is below the first grid point
+        } else if (upper_it == RT_TEMPERATURE_GRID_T9.end()) {
+            bound = BACK; // T9 is above the last grid point
+        } else {
+            bound = MIDDLE; // T9 is within the grid
+        }
+        const size_t upper_index = std::distance(RT_TEMPERATURE_GRID_T9.begin(), upper_it);
+        const size_t lower_index = upper_index - 1;
+        return {bound, data, upper_index, lower_index};
+    }
+
+    constexpr int RauscherThielemannPartitionFunction::make_key(
+        const int z,
+        const int a
+    ) {
+        return z * 1000 + a; // Simple key generation for Z and A
+    }
+}

src/network/meson.build

@@ -12,6 +12,9 @@ network_sources = files(
     'lib/screening/screening_types.cpp',
     'lib/screening/screening_weak.cpp',
     'lib/screening/screening_bare.cpp',
+    'lib/partition/partition_rauscher_thielemann.cpp',
+    'lib/partition/partition_ground.cpp',
+    'lib/partition/composite/partition_composite.cpp',
     'lib/utils/logging.cpp',
 )
 
@@ -58,6 +61,10 @@ network_headers = files(
     'include/gridfire/screening/screening_bare.h',
     'include/gridfire/screening/screening_weak.h',
     'include/gridfire/screening/screening_types.h',
+    'include/gridfire/partition/partition_abstract.h',
+    'include/gridfire/partition/partition_rauscher_thielemann.h',
+    'include/gridfire/partition/partition_ground.h',
+    'include/gridfire/partition/composite/partition_composite.h',
     'include/gridfire/utils/logging.h',
 )
 install_headers(network_headers, subdir : 'gridfire')

subprojects/libcomposition.wrap

@@ -1,4 +1,4 @@
 [wrap-git]
 url = https://github.com/4D-STAR/libcomposition.git
-revision = v1.1.0
+revision = v1.2.0
 depth = 1

tests/graphnet_sandbox/main.cpp

@@ -22,6 +22,8 @@
 #include <chrono>
 #include <functional>
 
+#include "gridfire/partition/composite/partition_composite.h"
+
 // Keep a copy of the previous handler
 static std::terminate_handler g_previousHandler = nullptr;
 
@@ -55,7 +57,7 @@ void quill_terminate_handler()
 int main() {
     g_previousHandler = std::set_terminate(quill_terminate_handler);
     quill::Logger* logger = fourdst::logging::LogManager::getInstance().getLogger("log");
-    logger->set_log_level(quill::LogLevel::Debug);
+    logger->set_log_level(quill::LogLevel::TraceL3);
     LOG_DEBUG(logger, "Starting Adaptive Engine View Example...");
 
     using namespace gridfire;
@@ -86,17 +88,25 @@ int main() {
     // }, "Approx8 Network Initialization");
     // std::cout << "Approx8 Network H-1: " << netOut.composition.getMassFraction("H-1") << " in " << netOut.num_steps << " steps." << std::endl;
 
+    using partition::BasePartitionType;
+    const auto partitionFunction = partition::CompositePartitionFunction({
+        BasePartitionType::RauscherThielemann,
+        BasePartitionType::GroundState
+    });
+    std::cout << "Partition Function for Mg-24: " << partitionFunction.evaluate(12, 24, 1.5) << std::endl;
+    std::cout << "Partition Function for F-23: " << partitionFunction.evaluate(9, 23, 1.5) << std::endl;
+    std::cout << "Partition Function for O-13: " << partitionFunction.evaluate(8, 13, 1.5) << std::endl;
 
     netIn.dt0 = 1e-15;
 
-    GraphEngine ReaclibEngine(composition);
-    ReaclibEngine.setPrecomputation(true);
-    // AdaptiveEngineView adaptiveEngine(ReaclibEngine);
-    io::SimpleReactionListFileParser parser{};
-    FileDefinedEngineView approx8EngineView(ReaclibEngine, "approx8.net", parser);
-    approx8EngineView.setScreeningModel(screening::ScreeningType::WEAK);
-    solver::QSENetworkSolver solver(approx8EngineView);
-    netOut = solver.evaluate(netIn);
+    // GraphEngine ReaclibEngine(composition);
+    // ReaclibEngine.setPrecomputation(true);
+    // // AdaptiveEngineView adaptiveEngine(ReaclibEngine);
+    // io::SimpleReactionListFileParser parser{};
+    // FileDefinedEngineView approx8EngineView(ReaclibEngine, "approx8.net", parser);
+    // approx8EngineView.setScreeningModel(screening::ScreeningType::WEAK);
+    // solver::QSENetworkSolver solver(approx8EngineView);
+    // netOut = solver.evaluate(netIn);
 
     // measure_execution_time([&]() {
     //     netOut = solver.evaluate(netIn);

utils/RauscherThielemann2000/format.py

@@ -0,0 +1,50 @@
+# generate_partition_binary.py
+import sys
+import struct
+import re
+
+def generate_binary_file(input_filepath, output_filepath):
+    RECORD_FORMAT = '<i i d 24d'
+
+    record_regex = re.compile(r"""
+        ^\s*[a-zA-Z0-9]+\s*\n
+        ^\s*
+        (?P<z>\d+)\s+
+        (?P<a>\d+)\s+
+        (?P<spin>[\d\.]+)
+        \s*\n
+        (?P<coeffs>
+          (?:^\s*(?:[\d.]+\s+){7}[\d.]+)\s*\n
+          (?:^\s*(?:[\d.]+\s+){7}[\d.]+)\s*\n
+          (?:^\s*(?:[\d.]+\s+){7}[\d.]+)
+        )
+    """, re.MULTILINE | re.VERBOSE)
+
+    with open(input_filepath, 'r') as f:
+        content = f.read()
+
+    record_count = 0
+    with open(output_filepath, 'wb') as f_out:
+        for match in record_regex.finditer(content):
+            z = int(match.group('z'))
+            a = int(match.group('a'))
+            spin = float(match.group('spin'))
+            coeffs_str = match.group('coeffs')
+            g_values = [float(val) for val in coeffs_str.split()]
+
+            if len(g_values) != 24:
+                print(f"Warning: Found {len(g_values)} coefficients for Z={z}, A={a}. Expected 24. Skipping.")
+                continue
+
+            packed_data = struct.pack(RECORD_FORMAT, z, a, spin, *g_values)
+            f_out.write(packed_data)
+            record_count += 1
+
+    print(f"Found and processed {record_count} records. Wrote to binary file {output_filepath}")
+
+
+if __name__ == "__main__":
+    if len(sys.argv) != 3:
+        print(f"Usage: python {sys.argv[0]} <input_partition_file> <output_binary_file>")
+        sys.exit(1)
+    generate_binary_file(sys.argv[1], sys.argv[2])