diff --git a/src/network/meson.build b/src/network/meson.build index 6b546f3f..ebf5448d 100644 --- a/src/network/meson.build +++ b/src/network/meson.build @@ -1,12 +1,14 @@ # Define the library network_sources = files( - 'private/network.cpp', - 'private/approx8.cpp' + 'private/network.cpp', + 'private/approx8.cpp', + 'private/netgraph.cpp' ) network_headers = files( 'public/network.h', - 'public/approx8.h' + 'public/approx8.h', + 'public/netgraph.h', ) dependencies = [ diff --git a/src/network/private/netgraph.cpp b/src/network/private/netgraph.cpp index ead0200b..adb84de5 100644 --- a/src/network/private/netgraph.cpp +++ b/src/network/private/netgraph.cpp @@ -1,3 +1,256 @@ -// -// Created by Emily Boudreaux on 6/19/25. -// +#include "netgraph.h" +#include "atomicSpecies.h" +#include "reaclib.h" +#include "network.h" +#include "species.h" + +#include "quill/LogMacros.h" + +#include +#include +#include +#include +#include +#include +#include +#include + + +namespace serif::network { + GraphNetwork::GraphNetwork( + const serif::composition::Composition &composition + ): + Network(REACLIB), + m_reactions(build_reaclib_nuclear_network(composition)), + m_initialComposition(composition), + m_currentComposition(composition), + m_cullingThreshold(0), + m_T9(0) { + syncInternalMaps(); + } + + GraphNetwork::GraphNetwork( + const serif::composition::Composition &composition, + const double cullingThreshold, + const double T9 + ): + Network(REACLIB), + m_reactions(build_reaclib_nuclear_network(composition, cullingThreshold, T9)), + m_initialComposition(composition), + m_currentComposition(composition), + m_cullingThreshold(cullingThreshold), + m_T9(T9) { + syncInternalMaps(); + } + + void GraphNetwork::syncInternalMaps() { + collectNetworkSpecies(); + populateReactionIDMap(); + populateSpeciesToIndexMap(); + } + + // --- Network Graph Construction Methods --- + void GraphNetwork::collectNetworkSpecies() { + m_networkSpecies.clear(); + m_networkSpeciesMap.clear(); + + std::set uniqueSpeciesNames; + + for (const auto& reaction: m_reactions) { + for (const auto& reactant: reaction.reactants()) { + uniqueSpeciesNames.insert(reactant.name()); + } + for (const auto& product: reaction.products()) { + uniqueSpeciesNames.insert(product.name()); + } + } + + for (const auto& name: uniqueSpeciesNames) { + auto it = serif::atomic::species.find(name); + if (it != serif::atomic::species.end()) { + m_networkSpecies.push_back(it->second); + m_networkSpeciesMap.insert({name, it->second}); + } else { + LOG_ERROR(m_logger, "Species '{}' not found in global atomic species database.", name); + throw std::runtime_error("Species not found in global atomic species database: " + std::string(name)); + } + } + + } + + void GraphNetwork::populateReactionIDMap() { + LOG_INFO(m_logger, "Populating reaction ID map for REACLIB graph network (serif::network::GraphNetwork)..."); + m_reactionIDMap.clear(); + for (const auto& reaction: m_reactions) { + m_reactionIDMap.insert({reaction.id(), reaction}); + } + LOG_INFO(m_logger, "Populated {} reactions in the reaction ID map.", m_reactionIDMap.size()); + } + + void GraphNetwork::populateSpeciesToIndexMap() { + m_speciesToIndexMap.clear(); + for (size_t i = 0; i < m_networkSpecies.size(); ++i) { + m_speciesToIndexMap.insert({m_networkSpecies[i], i}); + } + } + + void GraphNetwork::buildNetworkGraph() { + LOG_INFO(m_logger, "Building network graph..."); + m_reactions = build_reaclib_nuclear_network(m_initialComposition, m_cullingThreshold, m_T9); + syncInternalMaps(); + LOG_INFO(m_logger, "Network graph built with {} reactions and {} species.", m_reactions.size(), m_networkSpecies.size()); + } + // --- Basic Accessors and Queries --- + + const std::vector& GraphNetwork::getNetworkSpecies() const { + // Returns a constant reference to the vector of unique species in the network. + LOG_DEBUG(m_logger, "Providing access to network species vector. Size: {}.", m_networkSpecies.size()); + return m_networkSpecies; + } + + const reaclib::REACLIBReactionSet& GraphNetwork::getNetworkReactions() const { + // Returns a constant reference to the set of reactions in the network. + LOG_DEBUG(m_logger, "Providing access to network reactions set. Size: {}.", m_reactions.size()); + return m_reactions; + } + + bool GraphNetwork::involvesSpecies(const serif::atomic::Species& species) const { + // Checks if a given species is present in the network's species map for efficient lookup. + const bool found = m_networkSpeciesMap.contains(species.name()); + LOG_DEBUG(m_logger, "Checking if species '{}' is involved in the network: {}.", species.name(), found ? "Yes" : "No"); + return found; + } + + std::unordered_map GraphNetwork::getNetReactionStoichiometry(const reaclib::REACLIBReaction& reaction) const { + // Calculates the net stoichiometric coefficients for species in a given reaction. + std::unordered_map stoichiometry; + + // Iterate through reactants, decrementing their counts + for (const auto& reactant : reaction.reactants()) { + auto it = m_networkSpeciesMap.find(reactant.name()); + if (it != m_networkSpeciesMap.end()) { + stoichiometry[it->second]--; // Copy Species by value (PERF: Future performance improvements by using pointers or references) + } else { + LOG_WARNING(m_logger, "Reactant species '{}' in reaction '{}' not found in network species map during stoichiometry calculation.", + reactant.name(), reaction.id()); + } + } + + // Iterate through products, incrementing their counts + for (const auto& product : reaction.products()) { + auto it = m_networkSpeciesMap.find(product.name()); + if (it != m_networkSpeciesMap.end()) { + stoichiometry[it->second]++; // Copy Species by value (PERF: Future performance improvements by using pointers or references) + } else { + LOG_WARNING(m_logger, "Product species '{}' in reaction '{}' not found in network species map during stoichiometry calculation.", + product.name(), reaction.id()); + } + } + LOG_DEBUG(m_logger, "Calculated net stoichiometry for reaction '{}'. Total unique species in stoichiometry: {}.", reaction.id(), stoichiometry.size()); + return stoichiometry; + } + + // --- Validation Methods --- + bool GraphNetwork::validateConservation() const { + LOG_INFO(m_logger, "Validating mass (A) and charge (Z) conservation across all reactions in the network."); + + for (const auto& reaction : m_reactions) { + uint64_t totalReactantA = 0; + uint64_t totalReactantZ = 0; + uint64_t totalProductA = 0; + uint64_t totalProductZ = 0; + + // Calculate total A and Z for reactants + for (const auto& reactant : reaction.reactants()) { + auto it = m_networkSpeciesMap.find(reactant.name()); + if (it != m_networkSpeciesMap.end()) { + totalReactantA += it->second.a(); + totalReactantZ += it->second.z(); + } else { + // This scenario indicates a severe data integrity issue: + // a reactant is part of a reaction but not in the network's species map. + LOG_ERROR(m_logger, "CRITICAL ERROR: Reactant species '{}' in reaction '{}' not found in network species map during conservation validation.", + reactant.name(), reaction.id()); + return false; + } + } + + // Calculate total A and Z for products + for (const auto& product : reaction.products()) { + auto it = m_networkSpeciesMap.find(product.name()); + if (it != m_networkSpeciesMap.end()) { + totalProductA += it->second.a(); + totalProductZ += it->second.z(); + } else { + // Similar critical error for product species + LOG_ERROR(m_logger, "CRITICAL ERROR: Product species '{}' in reaction '{}' not found in network species map during conservation validation.", + product.name(), reaction.id()); + return false; + } + } + + // Compare totals for conservation + if (totalReactantA != totalProductA) { + LOG_ERROR(m_logger, "Mass number (A) not conserved for reaction '{}': Reactants A={} vs Products A={}.", + reaction.id(), totalReactantA, totalProductA); + return false; + } + if (totalReactantZ != totalProductZ) { + LOG_ERROR(m_logger, "Atomic number (Z) not conserved for reaction '{}': Reactants Z={} vs Products Z={}.", + reaction.id(), totalReactantZ, totalProductZ); + return false; + } + } + + LOG_INFO(m_logger, "Mass (A) and charge (Z) conservation validated successfully for all reactions."); + return true; // All reactions passed the conservation check + } + + // --- Generate Stoichiometry Matrix --- + void GraphNetwork::generateStoichiometryMatrix() { + LOG_INFO(m_logger, "Generating stoichiometry matrix..."); + + // Task 1: Set dimensions and initialize the matrix + size_t numSpecies = m_networkSpecies.size(); + size_t numReactions = m_reactions.size(); + m_stoichiometryMatrix.resize(numSpecies, numReactions, false); + + LOG_INFO(m_logger, "Stoichiometry matrix initialized with dimensions: {} rows (species) x {} columns (reactions).", + numSpecies, numReactions); + + // Task 2: Populate the stoichiometry matrix + // Iterate through all reactions, assign them a column index, and fill in their stoichiometric coefficients. + size_t reactionColumnIndex = 0; + for (const auto& reaction : m_reactions) { + // Get the net stoichiometry for the current reaction + std::unordered_map netStoichiometry = getNetReactionStoichiometry(reaction); + + // Iterate through the species and their coefficients in the stoichiometry map + for (const auto& pair : netStoichiometry) { + const serif::atomic::Species& species = pair.first; // The Species object + const int coefficient = pair.second; // The stoichiometric coefficient + + // Find the row index for this species + auto it = m_speciesToIndexMap.find(species); + if (it != m_speciesToIndexMap.end()) { + const size_t speciesRowIndex = it->second; + // Set the matrix element. Boost.uBLAS handles sparse insertion. + m_stoichiometryMatrix(speciesRowIndex, reactionColumnIndex) = coefficient; + } else { + // This scenario should ideally not happen if m_networkSpeciesMap and m_speciesToIndexMap are correctly synced + LOG_ERROR(m_logger, "CRITICAL ERROR: Species '{}' from reaction '{}' stoichiometry not found in species to index map.", + species.name(), reaction.id()); + throw std::runtime_error("Species not found in species to index map: " + std::string(species.name())); + } + } + reactionColumnIndex++; // Move to the next column for the next reaction + } + + LOG_INFO(m_logger, "Stoichiometry matrix population complete. Number of non-zero elements: {}.", + m_stoichiometryMatrix.nnz()); // Assuming nnz() exists for compressed_matrix + } + + NetOut GraphNetwork::evaluate(const NetIn &netIn) { + return Network::evaluate(netIn); + } +} diff --git a/src/network/public/netgraph.h b/src/network/public/netgraph.h index 0ceaff5d..759ad0b0 100644 --- a/src/network/public/netgraph.h +++ b/src/network/public/netgraph.h @@ -1,8 +1,63 @@ -// -// Created by Emily Boudreaux on 6/19/25. -// +#pragma once -#ifndef NETGRAPH_H -#define NETGRAPH_H +#include "network.h" +#include "reaclib.h" +#include "atomicSpecies.h" +#include "composition.h" -#endif //NETGRAPH_H +#include +#include +#include + +#include + + +// PERF: The function getNetReactionStoichiometry returns a map of species to their stoichiometric coefficients for a given reaction. +// this makes extra copies of the species, which is not ideal and could be optimized further. +// Even more relevant is the member m_reactionIDMap which makes copies of a REACLIBReaction for each reaction ID. +// REACLIBReactions are quite large data structures, so this could be a performance bottleneck. + +namespace serif::network { + class GraphNetwork final : public Network { + public: + explicit GraphNetwork(const serif::composition::Composition &composition); + explicit GraphNetwork(const serif::composition::Composition &composition, const double cullingThreshold, const double T9); + + NetOut evaluate(const NetIn &netIn) override; + + [[nodiscard]] const std::vector& getNetworkSpecies() const; + [[nodiscard]] const reaclib::REACLIBReactionSet& getNetworkReactions() const; + [[nodiscard]] std::unordered_map getNetReactionStoichiometry(const reaclib::REACLIBReaction& reaction) const; + + [[nodiscard]] bool involvesSpecies(const serif::atomic::Species& species) const; + + std::vector> detectCycles() const; + + std::vector topologicalSortReactions() const; + private: + reaclib::REACLIBReactionSet m_reactions; ///< Set of REACLIB reactions for the network. + serif::composition::Composition m_initialComposition; + serif::composition::Composition m_currentComposition; + double m_cullingThreshold; ///< Threshold for culling reactions. + double m_T9; ///< Temperature in T9 units to evaluate culling. + + std::vector m_networkSpecies; ///< The species in the network. + std::unordered_map m_networkSpeciesMap; + std::unordered_map m_reactionIDMap; ///< Map of reaction IDs to REACLIB reactions. + + boost::numeric::ublas::compressed_matrix m_stoichiometryMatrix; ///< Stoichiometry matrix for the network. + std::unordered_map m_speciesToIndexMap; ///< Map of species to their index in the stoichiometry matrix. + + + private: + void syncInternalMaps(); + void collectNetworkSpecies(); + void populateReactionIDMap(); + void populateSpeciesToIndexMap(); + void buildNetworkGraph(); + bool validateConservation() const; + void generateStoichiometryMatrix(); + + void pruneNetworkGraph(); + }; +}; \ No newline at end of file diff --git a/src/network/public/network.h b/src/network/public/network.h index 813e4888..4c9ca58f 100644 --- a/src/network/public/network.h +++ b/src/network/public/network.h @@ -33,11 +33,13 @@ namespace serif::network { enum NetworkFormat { APPROX8, ///< Approx8 nuclear reaction network format. + REACLIB, ///< General REACLIB nuclear reaction network format. UNKNOWN, }; static inline std::unordered_map FormatStringLookup = { {APPROX8, "Approx8"}, + {REACLIB, "REACLIB"}, {UNKNOWN, "Unknown"} };