Cache dependent

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()
@@ -870,6 +876,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

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
 

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@

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@")

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<r_prop_t>(); // reaction vertex property
+
+    std::set<v_desc_t> 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<stoic_t>(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<stoic_t>(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

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;
 };
 
 /**@}*/

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 <limits>
+#include <string>
+#include <iostream>
+#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<rrate_t>::max();
+  m_rmax = std::numeric_limits<rrate_t>::min();
+  m_rsum = static_cast<rrate_t>(0);
+  m_ravg = 0.0;
+  m_imin = std::numeric_limits<size_t>::max();
+  m_imax = std::numeric_limits<size_t>::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<double>(m_num_reactions);
+  m_iavg = m_isum / static_cast<double>(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__

src/reaction_network/reaction.hpp

@@ -29,6 +29,7 @@ class Reaction : public ReactionBase {
  public:
   using rdriver_t = std::template pair<VD, stoic_t>;
   using involved_species_t = std::template vector<rdriver_t>;
+  using vdesc_t = VD;
 
   Reaction();
   Reaction(const Reaction& rhs);
@@ -52,6 +53,12 @@ class Reaction : public ReactionBase {
     std::vector<std::string>& dep_params_nf);
 #endif // defined(WCS_HAS_EXPRTK)
 
+#ifdef WCS_CACHE_DEPENDENT
+  void set_dependent_reactions(const std::set<vdesc_t>& dependent);
+  const std::vector<vdesc_t>& 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<vdesc_t> m_dependent_reactions;
 };
 
 /**@}*/

src/reaction_network/reaction_impl.hpp

@@ -376,5 +376,25 @@ inline const typename Reaction<VD>::involved_species_t& Reaction<VD>::get_rate_i
   return m_rate_inputs;
 }
 
+#ifdef WCS_CACHE_DEPENDENT
+template <typename VD>
+void Reaction<VD>::set_dependent_reactions(const std::set<VD>& dependent)
+{
+  m_dependent_reactions.assign(dependent.begin(), dependent.end());
+}
+
+template <typename VD>
+const std::vector<VD>& Reaction<VD>::get_dependent_reactions() const
+{
+  return m_dependent_reactions;
+}
+
+template <typename VD>
+void Reaction<VD>::clear_dependent_reactions()
+{
+  m_dependent_reactions.clear();
+}
+#endif // WCS_CACHE_DEPENDENT
+
 /**@}*/
 } // end of namespace wcs

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