From 4ff0b6250cd866c9d66692c61020217656f49e40 Mon Sep 17 00:00:00 2001 From: Jae-Seung Yeom Date: Fri, 12 Feb 2021 03:03:47 -0800 Subject: [PATCH 1/4] added a function to cache dependent reactions --- src/reaction_network/network.cpp | 91 +++++++++++++++++++++- src/reaction_network/network.hpp | 10 +++ src/reaction_network/rate_stats.hpp | 101 +++++++++++++++++++++++++ src/reaction_network/reaction.hpp | 10 +++ src/reaction_network/reaction_impl.hpp | 20 +++++ 5 files changed, 231 insertions(+), 1 deletion(-) create mode 100644 src/reaction_network/rate_stats.hpp diff --git a/src/reaction_network/network.cpp b/src/reaction_network/network.cpp index 85a91a4..c50d374 100644 --- a/src/reaction_network/network.cpp +++ b/src/reaction_network/network.cpp @@ -190,6 +190,7 @@ void Network::init() m_reactions.reserve(num_vertices); m_species.reserve(num_vertices); + m_rstats.init(); v_iter_t vi, vi_end; for (boost::tie(vi, vi_end) = boost::vertices(m_graph); vi != vi_end; ++vi) { @@ -297,10 +298,12 @@ void Network::init() r.set_products(products); - set_reaction_rate(*vi); + // Set the reaction rate, and collect statistics. + m_rstats.read(set_reaction_rate(*vi), involved_species.size()); } } + m_rstats.set_avg(); sort_species(); build_index_maps(); @@ -825,5 +828,91 @@ void Network::print() const << num_inactive << "/" << reaction_list().size() << std::endl; } +#ifdef WCS_CACHE_DEPENDENT +void Network::cache_dependent_reactions( + std::function< bool(size_t, reaction_rate_t) >& to_cache) const +{ + const size_t n = m_my_reactions.size(); + + #if defined(_OPENMP) + #pragma omp parallel for //schedule(dynamic) + #endif // defined(_OPENMP) + for (size_t i = 0u; i < n; i++) { + const auto& rd = m_my_reactions[i]; // descriptor of reaction vertex + auto& rv = m_graph[rd]; // reaction vertex + auto& rp = rv.property(); // reaction vertex property + + std::set dependent_reactions; + + // ========================= reactant species ============================== + for (const auto ei_in : + boost::make_iterator_range(boost::in_edges(rd, m_graph))) + { + const auto sd_reac = boost::source(ei_in, m_graph); + #if !defined(__INTEL_COMPILER) + if constexpr (wcs::Vertex::_num_vertex_types_ > 3) { + // in case that there are other type of vertices than species or reaction + const auto& sv_reac = m_graph[sd_reac]; + if (sv_reac.get_type() != wcs::Vertex::_species_) continue; + } + #endif // !defined(__INTEL_COMPILER) + + const auto stoichio = m_graph[ei_in].get_stoichiometry_ratio(); + if (stoichio == static_cast(0)) { + continue; + } + + for (const auto vi_affected : + boost::make_iterator_range(boost::out_edges(sd_reac, m_graph))) + { + const auto rd_affected = boost::target(vi_affected, m_graph); + if (rd_affected == rd) continue; + if (m_graph[rd_affected].get_partition() != m_pid) continue; + dependent_reactions.insert(rd_affected); + } + } + + // ========================== product species ============================== + for (const auto ei_out : + boost::make_iterator_range(boost::out_edges(rd, m_graph))) + { + const auto sd_prod = boost::target(ei_out, m_graph); + #if !defined(__INTEL_COMPILER) + if constexpr (wcs::Vertex::_num_vertex_types_ > 3) { + // in case that there are other type of vertices than species or reaction + const auto& sv_prod = m_graph[sd_prod]; + if (sv_prod.get_type() != wcs::Vertex::_species_) continue; + } + #endif // !defined(__INTEL_COMPILER) + + const auto stoichio = m_graph[ei_out].get_stoichiometry_ratio(); + if (stoichio == static_cast(0)) { + continue; + } + + for (const auto vi_affected : + boost::make_iterator_range(boost::out_edges(sd_prod, m_graph))) + { + const auto rd_affected = boost::target(vi_affected, m_graph); + if (rd_affected == rd) continue; + if (m_graph[rd_affected].get_partition() != m_pid) continue; + dependent_reactions.insert(rd_affected); + } + } + + if (to_cache(dependent_reactions.size(), rp.get_rate())) { + rp.set_dependent_reactions(dependent_reactions); + } else { + rp.clear_dependent_reactions(); + } + } +} +#endif // WCS_CACHE_DEPENDENT + +const RateStats& Network::get_rate_stats() const +{ + return m_rstats; +} + /**@}*/ } // end of namespace wcs diff --git a/src/reaction_network/network.hpp b/src/reaction_network/network.hpp index e82a498..75e8e7b 100644 --- a/src/reaction_network/network.hpp +++ b/src/reaction_network/network.hpp @@ -32,6 +32,7 @@ #include "reaction_network/reaction.hpp" #include "reaction_network/vertex.hpp" #include "reaction_network/edge.hpp" +#include "reaction_network/rate_stats.hpp" namespace wcs { /** \addtogroup wcs_reaction_network @@ -166,6 +167,13 @@ class Network { /// Return the id of this partition partition_id_t get_partition_id() const; + #ifdef WCS_CACHE_DEPENDENT + void cache_dependent_reactions( + std::function< bool(size_t, reaction_rate_t) >& to_cache) const; + #endif // WCS_CACHE_DEPENDENT + + const RateStats& get_rate_stats() const; + void print() const; protected: @@ -220,6 +228,8 @@ class Network { */ rate_rules_dep_t m_rate_rules_dep_map; #endif // !defined(WCS_HAS_EXPRTK + + RateStats m_rstats; }; /**@}*/ diff --git a/src/reaction_network/rate_stats.hpp b/src/reaction_network/rate_stats.hpp new file mode 100644 index 0000000..376db21 --- /dev/null +++ b/src/reaction_network/rate_stats.hpp @@ -0,0 +1,101 @@ +/****************************************************************************** + * * + * Copyright 2020 Lawrence Livermore National Security, LLC and other * + * Whole Cell Simulator Project Developers. See the top-level COPYRIGHT * + * file for details. * + * * + * SPDX-License-Identifier: MIT * + * * + ******************************************************************************/ + +#ifndef __WCS_REACTION_NETWORK_RATE_STATS_HPP__ +#define __WCS_REACTION_NETWORK_RATE_STATS_HPP__ + +#include +#include +#include +#include "wcs_types.hpp" +#include "utils/to_string.hpp" + +namespace wcs { +/** \addtogroup wcs_reaction_network + * @{ */ + +struct RateStats { + using rrate_t = reaction_rate_t; + rrate_t m_rmin; + rrate_t m_rmax; + rrate_t m_rsum; + double m_ravg; + size_t m_imin; + size_t m_imax; + size_t m_isum; + double m_iavg; + size_t m_num_reactions; + + void init(); + void read(const rrate_t val, const size_t n_inputs); + void set_avg(); + void print() const; +}; + +inline void RateStats::init() +{ + m_rmin = std::numeric_limits::max(); + m_rmax = std::numeric_limits::min(); + m_rsum = static_cast(0); + m_ravg = 0.0; + m_imin = std::numeric_limits::max(); + m_imax = std::numeric_limits::min(); + m_isum = 0ul; + m_iavg = 0.0; + m_num_reactions = 0ul; +} + +inline void RateStats::read(const rrate_t val, const size_t num_inputs) +{ + if (val < m_rmin) m_rmin = val; + if (val > m_rmax) m_rmax = val; + m_rsum += val; + + if (num_inputs < m_imin) m_imin = num_inputs; + if (num_inputs > m_imax) m_imax = num_inputs; + m_isum += num_inputs; + + m_num_reactions ++; +} + +inline void RateStats::set_avg() +{ + m_ravg = m_rsum / static_cast(m_num_reactions); + m_iavg = m_isum / static_cast(m_num_reactions); +} + +inline void RateStats::print() const +{ + if (m_num_reactions == 0ul) { + return; + } + std::string msg = "Reaction statistic:"; + msg += "\n - min of rate: " + to_string_in_scientific(m_rmin); + msg += "\n - max of rate: " + to_string_in_scientific(m_rmax); + msg += "\n - sum of rate: " + to_string_in_scientific(m_rsum); + msg += "\n - avg of rate: " + to_string_in_scientific(m_ravg); + msg += "\n - min of number of inputs: " + std::to_string(m_imin); + msg += "\n - max of number of inputs: " + std::to_string(m_imax); + msg += "\n - sum of number of inputs: " + std::to_string(m_isum); + msg += "\n - avg of number of inputs: " + to_string_in_scientific(m_iavg); + msg += "\n - number of reactions: " + std::to_string(m_num_reactions); + msg += "\n"; + + #if defined(_OPENMP) + #pragma omp critical + #endif // defined(_OPENMP) + { + std::cout << msg; + } +} + +/**@}*/ +} // end of namespace wcs +#endif // __WCS_REACTION_NETWORK_RATE_STATSE_HPP__ diff --git a/src/reaction_network/reaction.hpp b/src/reaction_network/reaction.hpp index 53a4a3e..a71ed2c 100644 --- a/src/reaction_network/reaction.hpp +++ b/src/reaction_network/reaction.hpp @@ -29,6 +29,7 @@ class Reaction : public ReactionBase { public: using rdriver_t = std::template pair; using involved_species_t = std::template vector; + using vdesc_t = VD; Reaction(); Reaction(const Reaction& rhs); @@ -52,6 +53,12 @@ class Reaction : public ReactionBase { std::vector& dep_params_nf); #endif // defined(WCS_HAS_EXPRTK) +#ifdef WCS_CACHE_DEPENDENT + void set_dependent_reactions(const std::set& dependent); + const std::vector& get_dependent_reactions() const; + void clear_dependent_reactions(); +#endif // WCS_CACHE_DEPENDENT + protected: void reset(Reaction& obj); @@ -70,6 +77,9 @@ class Reaction : public ReactionBase { /// The BGL descriptors of input vertices to reaction rate formula involved_species_t m_rate_inputs; involved_species_t m_products; + + /// Cached list of dependent reactions that are local + std::vector m_dependent_reactions; }; /**@}*/ diff --git a/src/reaction_network/reaction_impl.hpp b/src/reaction_network/reaction_impl.hpp index 1033c4b..59efb75 100644 --- a/src/reaction_network/reaction_impl.hpp +++ b/src/reaction_network/reaction_impl.hpp @@ -376,5 +376,25 @@ inline const typename Reaction::involved_species_t& Reaction::get_rate_i return m_rate_inputs; } +#ifdef WCS_CACHE_DEPENDENT +template +void Reaction::set_dependent_reactions(const std::set& dependent) +{ + m_dependent_reactions.assign(dependent.begin(), dependent.end()); +} + +template +const std::vector& Reaction::get_dependent_reactions() const +{ + return m_dependent_reactions; +} + +template +void Reaction::clear_dependent_reactions() +{ + m_dependent_reactions.clear(); +} +#endif // WCS_CACHE_DEPENDENT + /**@}*/ } // end of namespace wcs From e1a65a94c22956f727350846500d412ae3c79430 Mon Sep 17 00:00:00 2001 From: Jae-Seung Yeom Date: Fri, 12 Feb 2021 03:06:31 -0800 Subject: [PATCH 2/4] CMake support for WCS_CACHE_DEPENDENT --- CMakeLists.txt | 1 + cmake/configure_files/WCSConfig.cmake.in | 1 + cmake/configure_files/wcs_config.hpp.in | 1 + cmake/configure_files/wcs_module.lua.in | 2 ++ src/ssa.cpp | 9 +++++++++ 5 files changed, 14 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 8e3c745..ebb5da0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -870,6 +870,7 @@ string(APPEND _str "\n") append_str_tf(_str WCS_GNU_LINUX WCS_64BIT_CNT + WCS_CACHE_DEPENDENT WCS_HAS_SUNDIALS WCS_HAS_SBML WCS_HAS_EXPRTK diff --git a/cmake/configure_files/WCSConfig.cmake.in b/cmake/configure_files/WCSConfig.cmake.in index d50eab6..4f2078b 100644 --- a/cmake/configure_files/WCSConfig.cmake.in +++ b/cmake/configure_files/WCSConfig.cmake.in @@ -50,6 +50,7 @@ set(WCS_HAS_METIS @WCS_HAS_METIS@) set(WCS_HAS_STD_FILESYSTEM @WCS_HAS_STD_FILESYSTEM@) set(WCS_HAS_PROTOBUF @WCS_HAS_PROTOBUF@) set(WCS_64BIT_CNT @WCS_64BIT_CNT@) +set(WCS_CACHE_DEPENDENT @WCS_CACHE_DEPENDENT@) # Setup dependencies diff --git a/cmake/configure_files/wcs_config.hpp.in b/cmake/configure_files/wcs_config.hpp.in index be44309..9dfc995 100644 --- a/cmake/configure_files/wcs_config.hpp.in +++ b/cmake/configure_files/wcs_config.hpp.in @@ -36,6 +36,7 @@ #cmakedefine WCS_HAS_STD_FILESYSTEM 1 #cmakedefine WCS_HAS_PROTOBUF 1 #cmakedefine WCS_64BIT_CNT 1 +#cmakedefine WCS_CACHE_DEPENDENT 1 #cmakedefine WCS_VERTEX_LIST_TYPE @WCS_VERTEX_LIST_TYPE@ #cmakedefine WCS_OUT_EDGE_LIST_TYPE @WCS_OUT_EDGE_LIST_TYPE@ diff --git a/cmake/configure_files/wcs_module.lua.in b/cmake/configure_files/wcs_module.lua.in index 452dcef..e81961b 100644 --- a/cmake/configure_files/wcs_module.lua.in +++ b/cmake/configure_files/wcs_module.lua.in @@ -24,6 +24,7 @@ -- WCS_HAS_STD_FILESYSTEM: @WCS_HAS_STD_FILESYSTEM@ -- WCS_HAS_PROTOBUF: @WCS_HAS_PROTOBUF@ -- WCS_64BIT_CNT: @WCS_64BIT_CNT@ +-- WCS_CACHE_DEPENDENT: @WCS_CACHE_DEPENDENT@ help( [[ @@ -58,6 +59,7 @@ whatis("WCS_HAS_METIS: @WCS_HAS_METIS@") whatis("WCS_HAS_STD_FILESYSTEM: @WCS_HAS_STD_FILESYSTEM@") whatis("WCS_HAS_PROTOBUF: @WCS_HAS_PROTOBUF@") whatis("WCS_64BIT_CNT: @WCS_64BIT_CNT@") +whatis("WCS_CACHE_DEPENDENT: @WCS_CACHE_DEPENDENT@") prepend_path("PATH","@CMAKE_INSTALL_PREFIX@/@CMAKE_INSTALL_BINDIR@") prepend_path("LD_LIBRARY_PATH","@CMAKE_INSTALL_PREFIX@/@CMAKE_INSTALL_LIBDIR@") diff --git a/src/ssa.cpp b/src/ssa.cpp index ae1de35..14c7958 100644 --- a/src/ssa.cpp +++ b/src/ssa.cpp @@ -37,6 +37,15 @@ int main(int argc, char** argv) wcs::Network& rnet = *rnet_ptr; rnet.load(cfg.m_infile); rnet.init(); + #ifdef WCS_CACHE_DEPENDENT + const auto& rstats = rnet.get_rate_stats(); + rstats.print(); + std::function< bool(size_t, wcs::reaction_rate_t) > to_cache + = [](size_t, wcs::reaction_rate_t) { + return true; + }; + rnet.cache_dependent_reactions(to_cache); + #endif // WCS_CACHE_DEPENDENT const wcs::Network::graph_t& g = rnet.graph(); if (!cfg.m_gvizfile.empty() && From 37f9bd1e00fee59af271a9f9f207bc81f89f83da Mon Sep 17 00:00:00 2001 From: Jae-Seung Yeom Date: Fri, 12 Feb 2021 09:50:17 -0800 Subject: [PATCH 3/4] Applied the dependent reaction caching to NRM --- src/sim_methods/CMakeLists.txt | 1 + src/sim_methods/fire_reaction_omp.hpp | 221 ++++++++++++++++++++++++++ src/sim_methods/sim_method.cpp | 111 ++++++++----- src/sim_methods/sim_state_change.hpp | 6 + src/sim_methods/ssa_nrm.cpp | 24 ++- src/sim_methods/ssa_nrm.hpp | 15 +- 6 files changed, 337 insertions(+), 41 deletions(-) create mode 100644 src/sim_methods/fire_reaction_omp.hpp diff --git a/src/sim_methods/CMakeLists.txt b/src/sim_methods/CMakeLists.txt index 618efa9..0cbc3c5 100644 --- a/src/sim_methods/CMakeLists.txt +++ b/src/sim_methods/CMakeLists.txt @@ -1,5 +1,6 @@ # Add the header and source files for this directory set_full_path(THIS_DIR_HEADERS + fire_reaction_omp.hpp sim_method.hpp sim_state_change.hpp ssa_nrm.hpp diff --git a/src/sim_methods/fire_reaction_omp.hpp b/src/sim_methods/fire_reaction_omp.hpp new file mode 100644 index 0000000..d404d87 --- /dev/null +++ b/src/sim_methods/fire_reaction_omp.hpp @@ -0,0 +1,221 @@ +/****************************************************************************** + * * + * Copyright 2020 Lawrence Livermore National Security, LLC and other * + * Whole Cell Simulator Project Developers. See the top-level COPYRIGHT * + * file for details. * + * * + * SPDX-License-Identifier: MIT * + * * + ******************************************************************************/ + +#ifndef WCS_SIM_METHODS_FIRE_REACTION_OMP_HPP +#define WCS_SIM_METHODS_FIRE_REACTION_OMP_HPP + +#if defined(_OPENMP) && \ + defined(WCS_CACHE_DEPENDENT) && \ + defined(WCS_OMP_REACTION_REACTANTS) && \ + defined(WCS_OMP_REACTION_PRODUCTS) && \ + !defined(ENABLE_SPECIES_UPDATE_TRACKING) + +#define WCS_OMP_FINE_GRAGIN 1 + +#include +#include +#include "sim_methods/sim_method.hpp" + + +namespace wcs { +/** \addtogroup wcs_reaction_network + * @{ */ + +template +void reduce_union(typename std::vector< std::vector >& results) +{ + const size_t num_srcs = results.size(); + size_t n = results.size(); + size_t d = 1ul; + + while ((n/=2)) + { + #pragma omp parallel for + for (size_t i = 0ul; i < num_srcs; i += 2*d) + { + if (i+d >= num_srcs) { + continue; + } + const typename std::vector& src1 = results[i]; + const typename std::vector& src2 = results[i+d]; + typename std::vector uni(src1.size() + src2.size()); + typename std::vector::iterator it; + + it = std::set_union(src1.begin(), src1.end(), + src2.begin(), src2.end(), + uni.begin()); + uni.resize(it-uni.begin()); + results[i] = std::move(uni); + results[i+d].clear(); + } + d *= 2; + } + + if (num_srcs % 2) { + const typename std::vector& src1 = results[0]; + const typename std::vector& src2 = results[d]; + typename std::vector uni(src1.size() + src2.size()); + typename std::vector::iterator it; + + it = std::set_union(src1.begin(), src1.end(), + src2.begin(), src2.end(), + uni.begin()); + uni.resize(it-uni.begin()); + results[0].swap(uni); + } +} + +/** + * Execute the chosen reaction. + * In other words, update the species population involved in the reaction. + * In addition, record how the species are updated, and which other reactions + * are affected as a result. + * The former can be used to undo the reaction if needed. The latter is used + * to update the propensity of the reactions affected by the changes in species + * counts. + */ +bool Sim_Method::fire_reaction(Sim_State_Change& digest) +{ + using s_prop_t = wcs::Species; + using v_desc_t = wcs::Network::v_desc_t; + + // Reactions do not change the connectivity, but only change the property + // data, which are allocated outside of BGL graph, but only linked to it. + const wcs::Network::graph_t& g = m_net_ptr->graph(); + // The vertex descriptor of the reaction to undo + const auto& rd_firing = digest.m_reaction_fired; + const auto& rv_firing = g[rd_firing]; + const auto& rp_firing = rv_firing.property(); + const bool not_cached = rp_firing.get_dependent_reactions().empty(); + std::vector< std::vector > reactions_affected(m_num_threads); + + + #if defined(WCS_OMP_RUN_PARTITION) + const auto pid = m_net_ptr->get_partition_id(); + #endif // defined(WCS_OMP_RUN_PARTITION) + + // ========================= reactant species ================================ + using e_desc_t = wcs::Network::e_desc_t; + using ie_iter_t = boost::graph_traits::in_edge_iterator; + ie_iter_t iei, iei_end; + + boost::tie(iei, iei_end) = boost::in_edges(rd_firing, g); + const std::vector inedges(iei, iei_end); + const size_t nie = inedges.size(); + + #pragma omp parallel for //schedule(dynamic) + for (size_t i = 0u; i < nie; ++i) + { + const auto& ei_in = inedges[i]; + const auto vd_updating = boost::source(ei_in, g); + const auto& sv_updating = g[vd_updating]; + if constexpr (wcs::Vertex::_num_vertex_types_ > 3) { + // in case that there are other type of vertices than species or reaction + if (sv_updating.get_type() != wcs::Vertex::_species_) continue; + } + + auto& sp_updating = sv_updating.property(); + const auto stoichio = g[ei_in].get_stoichiometry_ratio(); + if (stoichio == static_cast(0)) { + continue; + } + #ifdef NDEBUG + sp_updating.dec_count(stoichio); + #else + // This really should not happen because whether the reaction is feasible is + // checked before computing reaction time or propensity. + if (!sp_updating.dec_count(stoichio)) { // State update + std::string err = "Not enough reactants of " + sv_updating.get_label() + + "[" + std::to_string(sp_updating.get_count()) + + "] for reaction " + g[rd_firing].get_label(); + WCS_THROW(err); + //return false; // TODO: graceful termination + } + #endif + + if (not_cached) + { + for (const auto vi_affected : + boost::make_iterator_range(boost::out_edges(vd_updating, g))) + { + const auto rd_affected = boost::target(vi_affected, g); + if (rd_affected == rd_firing) continue; + #if defined(WCS_OMP_RUN_PARTITION) + if (m_net_ptr->graph()[rd_affected].get_partition() != pid) continue; + #endif // defined(WCS_OMP_RUN_PARTITION) + reactions_affected[omp_get_thread_num()].push_back(rd_affected); + } + } + } + + // ========================== product species ================================ + using oe_iter_t = boost::graph_traits::out_edge_iterator; + oe_iter_t oei, oei_end; + boost::tie(oei, oei_end) = boost::out_edges(rd_firing, g); + const std::vector outedges(oei, oei_end); + const size_t noe = outedges.size(); + + #pragma omp parallel for //schedule(dynamic) + for (size_t i = 0u; i < noe; ++i) + { + const auto& ei_out = outedges[i]; + const auto vd_updating = boost::target(ei_out, g); + const auto& sv_updating = g[vd_updating]; + if constexpr (wcs::Vertex::_num_vertex_types_ > 3) { + // in case that there are other type of vertices than species or reaction + if (sv_updating.get_type() != wcs::Vertex::_species_) continue; + } + + auto& sp_updating = sv_updating.property(); + const auto stoichio = g[ei_out].get_stoichiometry_ratio(); + if (stoichio == static_cast(0)) { + continue; + } + #ifdef NDEBUG + sp_updating.inc_count(stoichio); + #else + if (!sp_updating.inc_count(stoichio)) { // State update + std::string err = "Can not produce more of " + sv_updating.get_label() + + "[" + std::to_string(sp_updating.get_count()) + + "] by reaction " + g[rd_firing].get_label(); + WCS_THROW(err); + //return false; // TODO: graceful termination + } + #endif + + if (not_cached) + { + for (const auto vi_affected : + boost::make_iterator_range(boost::out_edges(vd_updating, g))) + { + const auto rd_affected = boost::target(vi_affected, g); + if (rd_affected == rd_firing) continue; + #if defined(WCS_OMP_RUN_PARTITION) + if (m_net_ptr->graph()[rd_affected].get_partition() != pid) continue; + #endif // defined(_OPENMP) && defined(WCS_OMP_RUN_PARTITION) + reactions_affected[omp_get_thread_num()].push_back(rd_affected); + } + } + } + + if (not_cached) { + reduce_union(reactions_affected); + digest.m_dependent_reactions.assign(reactions_affected[0].begin(), reactions_affected[0].end()); + } else { + digest.m_dependent_reactions = rp_firing.get_dependent_reactions(); + } + return true; +} + +/**@}*/ +} // end of namespace wcs + +#endif //defined(_OPENMP) +#endif // WCS_SIM_METHODS_FIRE_REACTION_OMP_HPP diff --git a/src/sim_methods/sim_method.cpp b/src/sim_methods/sim_method.cpp index 6eca129..f2c7c8a 100644 --- a/src/sim_methods/sim_method.cpp +++ b/src/sim_methods/sim_method.cpp @@ -11,6 +11,7 @@ #include #include // std::forward #include "sim_methods/sim_method.hpp" +#include "sim_methods/fire_reaction_omp.hpp" namespace wcs { /** \addtogroup wcs_reaction_network @@ -110,6 +111,7 @@ void Sim_Method::finalize_recording() { } } +#if !WCS_OMP_FINE_GRAGIN /** * Execute the chosen reaction. * In other words, update the species population involved in the reaction. @@ -128,14 +130,21 @@ bool Sim_Method::fire_reaction(Sim_State_Change& digest) const wcs::Network::graph_t& g = m_net_ptr->graph(); // The vertex descriptor of the reaction to undo const auto& rd_firing = digest.m_reaction_fired; + #ifdef WCS_CACHE_DEPENDENT + const auto& rv_firing = g[rd_firing]; + const auto& rp_firing = rv_firing.property(); + const bool not_cached = rp_firing.get_dependent_reactions().empty(); + std::set reactions_affected; + #else + auto& reactions_affected = digest.m_reactions_affected; + reactions_affected.clear(); + #endif // WCS_CACHE_DEPENDENT + #ifdef ENABLE_SPECIES_UPDATE_TRACKING // can be used for undo_species_updates() auto& updating_species = digest.m_species_updated; updating_species.clear(); #endif // ENABLE_SPECIES_UPDATE_TRACKING - auto& reactions_affected = digest.m_reactions_affected; - - reactions_affected.clear(); #if defined(_OPENMP) && defined(WCS_OMP_RUN_PARTITION) const auto pid = m_net_ptr->get_partition_id(); @@ -187,17 +196,22 @@ bool Sim_Method::fire_reaction(Sim_State_Change& digest) } #endif // ENABLE_SPECIES_UPDATE_TRACKING - for (const auto vi_affected : - boost::make_iterator_range(boost::out_edges(vd_updating, g))) + #ifdef WCS_CACHE_DEPENDENT + if (not_cached) + #endif // WCS_CACHE_DEPENDENT { - const auto rd_affected = boost::target(vi_affected, g); - if (rd_affected == rd_firing) continue; - #if defined(_OPENMP) && defined(WCS_OMP_RUN_PARTITION) - if (m_net_ptr->graph()[rd_affected].get_partition() != pid) continue; - #endif // defined(_OPENMP) && defined(WCS_OMP_RUN_PARTITION) - #pragma omp critical + for (const auto vi_affected : + boost::make_iterator_range(boost::out_edges(vd_updating, g))) { - reactions_affected.insert(rd_affected); + const auto rd_affected = boost::target(vi_affected, g); + if (rd_affected == rd_firing) continue; + #if defined(_OPENMP) && defined(WCS_OMP_RUN_PARTITION) + if (m_net_ptr->graph()[rd_affected].get_partition() != pid) continue; + #endif // defined(_OPENMP) && defined(WCS_OMP_RUN_PARTITION) + #pragma omp critical + { + reactions_affected.insert(rd_affected); + } } } } @@ -235,15 +249,20 @@ bool Sim_Method::fire_reaction(Sim_State_Change& digest) updating_species.emplace_back(std::make_pair(vd_updating, -stoichio)); #endif // ENABLE_SPECIES_UPDATE_TRACKING - for (const auto vi_affected : - boost::make_iterator_range(boost::out_edges(vd_updating, g))) + #ifdef WCS_CACHE_DEPENDENT + if (not_cached) + #endif // WCS_CACHE_DEPENDENT { - const auto rd_affected = boost::target(vi_affected, g); - if (rd_affected == rd_firing) continue; - #if defined(_OPENMP) && defined(WCS_OMP_RUN_PARTITION) - if (m_net_ptr->graph()[rd_affected].get_partition() != pid) continue; - #endif // defined(_OPENMP) && defined(WCS_OMP_RUN_PARTITION) - reactions_affected.insert(rd_affected); + for (const auto vi_affected : + boost::make_iterator_range(boost::out_edges(vd_updating, g))) + { + const auto rd_affected = boost::target(vi_affected, g); + if (rd_affected == rd_firing) continue; + #if defined(_OPENMP) && defined(WCS_OMP_RUN_PARTITION) + if (m_net_ptr->graph()[rd_affected].get_partition() != pid) continue; + #endif // defined(_OPENMP) && defined(WCS_OMP_RUN_PARTITION) + reactions_affected.insert(rd_affected); + } } } #endif // defined(_OPENMP) && defined(WCS_OMP_REACTION_REACTANTS) // ---------- @@ -290,17 +309,22 @@ bool Sim_Method::fire_reaction(Sim_State_Change& digest) } #endif // ENABLE_SPECIES_UPDATE_TRACKING - for (const auto vi_affected : - boost::make_iterator_range(boost::out_edges(vd_updating, g))) + #ifdef WCS_CACHE_DEPENDENT + if (not_cached) + #endif // WCS_CACHE_DEPENDENT { - const auto rd_affected = boost::target(vi_affected, g); - if (rd_affected == rd_firing) continue; - #if defined(_OPENMP) && defined(WCS_OMP_RUN_PARTITION) - if (m_net_ptr->graph()[rd_affected].get_partition() != pid) continue; - #endif // defined(_OPENMP) && defined(WCS_OMP_RUN_PARTITION) - #pragma omp critical + for (const auto vi_affected : + boost::make_iterator_range(boost::out_edges(vd_updating, g))) { - reactions_affected.insert(rd_affected); + const auto rd_affected = boost::target(vi_affected, g); + if (rd_affected == rd_firing) continue; + #if defined(_OPENMP) && defined(WCS_OMP_RUN_PARTITION) + if (m_net_ptr->graph()[rd_affected].get_partition() != pid) continue; + #endif // defined(_OPENMP) && defined(WCS_OMP_RUN_PARTITION) + #pragma omp critical + { + reactions_affected.insert(rd_affected); + } } } } @@ -337,21 +361,34 @@ bool Sim_Method::fire_reaction(Sim_State_Change& digest) updating_species.emplace_back(std::make_pair(vd_updating, stoichio)); #endif // ENABLE_SPECIES_UPDATE_TRACKING - for (const auto vi_affected : - boost::make_iterator_range(boost::out_edges(vd_updating, g))) + #ifdef WCS_CACHE_DEPENDENT + if (not_cached) + #endif // WCS_CACHE_DEPENDENT { - const auto rd_affected = boost::target(vi_affected, g); - if (rd_affected == rd_firing) continue; - #if defined(_OPENMP) && defined(WCS_OMP_RUN_PARTITION) - if (m_net_ptr->graph()[rd_affected].get_partition() != pid) continue; - #endif // defined(_OPENMP) && defined(WCS_OMP_RUN_PARTITION) - reactions_affected.insert(rd_affected); + for (const auto vi_affected : + boost::make_iterator_range(boost::out_edges(vd_updating, g))) + { + const auto rd_affected = boost::target(vi_affected, g); + if (rd_affected == rd_firing) continue; + #if defined(_OPENMP) && defined(WCS_OMP_RUN_PARTITION) + if (m_net_ptr->graph()[rd_affected].get_partition() != pid) continue; + #endif // defined(_OPENMP) && defined(WCS_OMP_RUN_PARTITION) + reactions_affected.insert(rd_affected); + } } } #endif // defined(_OPENMP) && defined(WCS_OMP_REACTION_PRODUCTS) // ----------- + #ifdef WCS_CACHE_DEPENDENT + if (!not_cached) { + digest.m_dependent_reactions = rp_firing.get_dependent_reactions(); + } else { + digest.m_dependent_reactions.assign(reactions_affected.begin(), reactions_affected.end()); + } + #endif // WCS_CACHE_DEPENDENT return true; } +#endif // !WCS_OMP_FINE_GRAGIN #ifdef ENABLE_SPECIES_UPDATE_TRACKING diff --git a/src/sim_methods/sim_state_change.hpp b/src/sim_methods/sim_state_change.hpp index 16b675e..f475947 100644 --- a/src/sim_methods/sim_state_change.hpp +++ b/src/sim_methods/sim_state_change.hpp @@ -29,6 +29,11 @@ struct Sim_State_Change { /** Type for the list of reactions that share any of the species with the * firing reaction */ using affected_reactions_t = std::set; + /** To carry the same information as affected_reactions_t does. However, + * by using the randomly accessible container, this makes it easier to + * construct a parallel loop over the container. + */ + using dependent_reactions_t = std::vector; using reaction_times_t = std::vector >; using revent_t = std::pair; @@ -48,6 +53,7 @@ struct Sim_State_Change { * result of processing the reaction. */ affected_reactions_t m_reactions_affected; + dependent_reactions_t m_dependent_reactions; /** * The currently expected times of the affected reactions to occur. diff --git a/src/sim_methods/ssa_nrm.cpp b/src/sim_methods/ssa_nrm.cpp index 77c590e..c509bdf 100644 --- a/src/sim_methods/ssa_nrm.cpp +++ b/src/sim_methods/ssa_nrm.cpp @@ -229,7 +229,11 @@ wcs::sim_time_t SSA_NRM::adjust_reaction_time(const v_desc_t& vd, * this function is called when the reaction fired is local. */ void SSA_NRM::update_reactions(const sim_time_t t_fired, - const Sim_Method::affected_reactions_t& affected, + #ifdef WCS_CACHE_DEPENDENT + const dependent_reactions_t& affected, + #else + const affected_reactions_t& affected, + #endif // WCS_CACHE_DEPENDENT SSA_NRM::reaction_times_t& affected_rtimes) { #if defined(_OPENMP) && defined(WCS_OMP_RUN_PARTITION) @@ -239,7 +243,11 @@ void SSA_NRM::update_reactions(const sim_time_t t_fired, lambdas_for_indexed_heap #if defined(_OPENMP) && defined(WCS_OMP_REACTION_UPDATES) + #ifdef WCS_CACHE_DEPENDENT + const std::vector& r_affected = affected; + #else const std::vector r_affected(affected.begin(), affected.end()); + #endif // WCS_CACHE_DEPENDENT #pragma omp parallel for for (size_t i = 0ul; i < r_affected.size(); i++) { @@ -277,7 +285,11 @@ void SSA_NRM::update_reactions(const sim_time_t t_fired, */ void SSA_NRM::update_reactions( const SSA_NRM::priority_t& fired, - const Sim_Method::affected_reactions_t& affected, + #ifdef WCS_CACHE_DEPENDENT + const dependent_reactions_t& affected, + #else + const affected_reactions_t& affected, + #endif // WCS_CACHE_DEPENDENT SSA_NRM::reaction_times_t& affected_rtimes) { const auto& t_fired = fired.first; @@ -299,7 +311,11 @@ void SSA_NRM::update_reactions( #endif // defined(WCS_HAS_ROSS) #if defined(_OPENMP) && defined(WCS_OMP_REACTION_UPDATES) + #ifdef WCS_CACHE_DEPENDENT + const std::vector& r_affected = affected; + #else const std::vector r_affected(affected.begin(), affected.end()); + #endif // WCS_CACHE_DEPENDENT #pragma omp parallel for for (size_t i = 0ul; i < r_affected.size(); i++) { @@ -484,7 +500,11 @@ bool SSA_NRM::forward(const revent_t firing) Sim_Method::fire_reaction(digest); // update the propensities and times of those reactions fired and affected + #ifdef WCS_CACHE_DEPENDENT + update_reactions(firing, digest.m_dependent_reactions, digest.m_reaction_times); + #else update_reactions(firing, digest.m_reactions_affected, digest.m_reaction_times); + #endif // WCS_CACHE_DEPENDENT #if !defined(WCS_HAS_ROSS) // With ROSS, tracing and sampling are moved to process at commit time diff --git a/src/sim_methods/ssa_nrm.hpp b/src/sim_methods/ssa_nrm.hpp index f9419eb..b94d411 100644 --- a/src/sim_methods/ssa_nrm.hpp +++ b/src/sim_methods/ssa_nrm.hpp @@ -36,6 +36,9 @@ class SSA_NRM : public Sim_Method { /// Type of the pair of BGL vertex descriptor for reaction and the its time using reaction_times_t = Sim_State_Change::reaction_times_t; + using affected_reactions_t = Sim_State_Change::affected_reactions_t; + using dependent_reactions_t = Sim_State_Change::dependent_reactions_t; + SSA_NRM(const std::shared_ptr& net_ptr); SSA_NRM(SSA_NRM&& other) = default; @@ -83,12 +86,20 @@ class SSA_NRM : public Sim_Method { /// Used in parallel mode with partitioned network bool advance_time_and_iter(const sim_time_t t_new); void update_reactions(const sim_time_t t_fired, - const Sim_Method::affected_reactions_t& affected, + #ifdef WCS_CACHE_DEPENDENT + const dependent_reactions_t& affected, + #else + const affected_reactions_t& affected, + #endif // WCS_CACHE_DEPENDENT reaction_times_t& affected_rtimes); #endif // defined(_OPENMP) && defined(WCS_OMP_RUN_PARTITION) void update_reactions(const priority_t& fired, - const Sim_Method::affected_reactions_t& affected, + #ifdef WCS_CACHE_DEPENDENT + const dependent_reactions_t& affected, + #else + const affected_reactions_t& affected, + #endif // WCS_CACHE_DEPENDENT reaction_times_t& affected_rtimes); protected: From 1accd0d7a5bf5cdc28ece8afc6c9de2188791afd Mon Sep 17 00:00:00 2001 From: Jae-Seung Yeom Date: Fri, 12 Feb 2021 09:51:14 -0800 Subject: [PATCH 4/4] use it in the SSA driver --- CMakeLists.txt | 6 ++++++ src/ssa_partition.cpp | 18 ++++++++++++++++++ 2 files changed, 24 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index ebb5da0..7762109 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -39,6 +39,7 @@ string(TOLOWER "${PROJECT_NAME}" LOWER_PROJECT_NAME) set(ignoreUnused ${WCS_OMP_MODE}${WCS_WITH_VTUN}${VTUNE_DIR}) set(ignoreUnused ${WCS_WITH_HPCTOOLKIT}${HPCTOOLKIT_DIR}) set(ignoreUnused ${WCS_WITH_ROSS}${ROSS_ROOT}${ROSS_DIR}) +set(ignoreUnused ${WCS_CACHE_DEPENDENT}) ################################################################ # Version setup @@ -373,6 +374,11 @@ if (OpenMP_CXX_FOUND) elseif(WCS_OMP_MODE MATCHES 5) set(WCS_OMP_REACTION_UPDATES ON) set(WCS_OMP_RUN_PARTITION ON) + elseif(WCS_OMP_MODE MATCHES 6) + set(WCS_OMP_REACTION_UPDATES ON) + set(WCS_OMP_REACTION_REACTANTS ON) + set(WCS_OMP_REACTION_PRODUCTS ON) + set(WCS_OMP_RUN_PARTITION ON) elseif() message(STATUS "Unknown WCS_OMP_MODE (${WCS_OMP_MODE})") endif() diff --git a/src/ssa_partition.cpp b/src/ssa_partition.cpp index 5830906..be9f450 100644 --- a/src/ssa_partition.cpp +++ b/src/ssa_partition.cpp @@ -250,6 +250,16 @@ void wcs_init(const wcs::SSA_Params& cfg, const partition_idx_t& parts) net.init(); net.set_partition(parts, tid); + #ifdef WCS_CACHE_DEPENDENT + const auto& rstats = net.get_rate_stats(); + rstats.print(); + std::function< bool(size_t, wcs::reaction_rate_t) > to_cache + = [](size_t, wcs::reaction_rate_t) { + return true; + }; + net.cache_dependent_reactions(to_cache); + #endif // WCS_CACHE_DEPENDENT + ssa_ptr = std::make_unique(net_ptr); wcs::SSA_NRM& ssa = *ssa_ptr; ssa.set_num_threads(shared_state.m_num_inner_threads); @@ -395,10 +405,18 @@ bool forward(const nrm_evt_t& evt_earliest) ssa.fire_reaction(digest); if (local) { // Update the propensities and times of all local reactions that are fired and affected + #ifdef WCS_CACHE_DEPENDENT + ssa.update_reactions(firing, digest.m_dependent_reactions, digest.m_reaction_times); + #else ssa.update_reactions(firing, digest.m_reactions_affected, digest.m_reaction_times); + #endif // WCS_CACHE_DEPENDENT } else { // This does not update the reaction fired which is not local. + #ifdef WCS_CACHE_DEPENDENT + ssa.update_reactions(firing.first, digest.m_dependent_reactions, digest.m_reaction_times); + #else ssa.update_reactions(firing.first, digest.m_reactions_affected, digest.m_reaction_times); + #endif // WCS_CACHE_DEPENDENT } #pragma omp master