From 1aea1c4e11a8baa73bb1b193f06429001d12e71e Mon Sep 17 00:00:00 2001 From: Konstantia Georgouli Date: Thu, 27 May 2021 20:26:13 -0700 Subject: [PATCH 01/14] Support initialConcentration, constant reactants and units --- src/reaction_network/network.cpp | 70 ++++++-- src/reaction_network/vertex.hpp | 111 +++++++++--- src/utils/generate_cxx_code.cpp | 280 +++++++++++++++++++++++++++++-- src/utils/graph_factory.hpp | 78 ++++++--- src/utils/sbml_utils.cpp | 30 +++- 5 files changed, 483 insertions(+), 86 deletions(-) diff --git a/src/reaction_network/network.cpp b/src/reaction_network/network.cpp index 85a91a4..43eca39 100644 --- a/src/reaction_network/network.cpp +++ b/src/reaction_network/network.cpp @@ -22,6 +22,11 @@ #include // is_same<> #include // lexicographical_compare(), sort() #include // numeric_limits +#include + +#define TIME_DIFF(Tstart, Tend ) \ + ((double) (1000000000L * ((Tend).tv_sec - (Tstart).tv_sec) + \ + (Tend).tv_nsec - (Tstart).tv_nsec) / 1000000000L) #if defined(WCS_HAS_SBML) #include @@ -154,9 +159,13 @@ void Network::loadSBML(const std::string sbml_filename, const bool reuse) generate_cxx_code code_generator(lib_filename, !reuse); typename params_map_t::const_iterator pit; + struct timespec t_1, t_2; + clock_gettime(CLOCK_MONOTONIC_RAW, &t_1); code_generator.generate_code(*model, m_dep_params_f, m_dep_params_nf, m_rate_rules_dep_map); - + clock_gettime(CLOCK_MONOTONIC_RAW, &t_2); + using std::operator<<; + std::cout << "Time for creating the generated file: " << TIME_DIFF(t_1, t_2 ) << " sec" << std::endl; const std::string library_file = code_generator.compile_code(); using std::operator<<; @@ -192,12 +201,17 @@ void Network::init() m_species.reserve(num_vertices); v_iter_t vi, vi_end; + int reaction_sz = 0; + double sum_in = 0, sum_reac_mod = 0, sum_involved = 0, sum_out=0; + for (boost::tie(vi, vi_end) = boost::vertices(m_graph); vi != vi_end; ++vi) { const v_prop_t& v = m_graph[*vi]; const auto vt = static_cast(v.get_typeid()); if (vt == v_prop_t::_species_) { m_species.emplace_back(*vi); } else { + reaction_sz = reaction_sz +1; + int reaction_in = 0, reac_mod_in = 0, reaction_out =0; using directed_category = boost::graph_traits::directed_category; constexpr bool is_bidirectional = std::is_same::value; @@ -227,17 +241,22 @@ void Network::init() boost::make_iterator_range(boost::in_edges(reaction, m_graph))) { v_desc_t reactant = boost::source(ei_in, m_graph); - #if !defined(WCS_HAS_EXPRTK) - //check for the reactants which are not actually reactants and put their stoichiometry 0 - if ( std::find(params_reactants.begin(), params_reactants.end(), - m_graph[reactant].get_label()) == params_reactants.end()) { - //std::cout << "Not reactant " << m_graph[reactant].get_label() << std::endl; - m_graph[ei_in].set_stoichiometry_ratio(0); - } - #endif // !defined(WCS_HAS_EXPRTK) + // #if !defined(WCS_HAS_EXPRTK) + // //check for the reactants which are not actually reactants and put their stoichiometry 0 + // if ( std::find(params_reactants.begin(), params_reactants.end(), + // m_graph[reactant].get_label()) == params_reactants.end()) { + // //std::cout << "Not reactant " << m_graph[reactant].get_label() << std::endl; + // m_graph[ei_in].set_stoichiometry_ratio(0); + // } + // #endif // !defined(WCS_HAS_EXPRTK) const auto st = m_graph[ei_in].get_stoichiometry_ratio(); involved_species.insert(std::make_pair(m_graph[reactant].get_label(), std::make_pair(reactant, st))); + + if (m_graph[ei_in].get_stoichiometry_ratio() > 0 ){ + reaction_in = reaction_in + 1; + } + reac_mod_in = reac_mod_in + 1; } } @@ -246,6 +265,7 @@ void Network::init() v_desc_t product = boost::target(ei_out, m_graph); products.insert(std::make_pair(m_graph[product].get_label(), std::make_pair(product, 1))); + reaction_out = reaction_out + 1; } m_reactions.emplace_back(reaction); @@ -295,12 +315,31 @@ void Network::init() r.set_rate_inputs(involved_species); #endif // !defined(WCS_HAS_EXPRTK) + auto & rprop = m_graph[*vi].checked_property< Reaction >(); + const auto& ri = rprop.get_rate_inputs(); + + sum_involved = sum_involved + ri.size(); + sum_in = sum_in + reaction_in; + sum_reac_mod = sum_reac_mod + reac_mod_in; + sum_out = sum_out + reaction_out; + //using std::operator<<; + //std::cout << "involved species: " << ri.size() << " reactants:" << reaction_in << std::endl; + r.set_products(products); set_reaction_rate(*vi); } } + double mean_in = sum_in / reaction_sz; + double mean_out = sum_out / reaction_sz; + double mean_reac_mod = sum_reac_mod / reaction_sz; + double mean_involved = sum_involved / reaction_sz; + using std::operator<<; + + std::cout << "Reactions " << reaction_sz << " mean of inputs " << mean_reac_mod + << ", mean of reactants " << mean_in << ", mean of species involved in reaction rate " << mean_involved + << " and, mean of products " << mean_out << std::endl; sort_species(); build_index_maps(); @@ -351,7 +390,9 @@ reaction_rate_t Network::set_reaction_rate(const Network::v_desc_t r) const double Network::compute_all_reaction_rates(const unsigned n) const { - double t_start = get_time(); + //double t_start = get_time(); + struct timespec t_start, t_end; + clock_gettime(CLOCK_MONOTONIC_RAW, &t_start); for (unsigned i = 0u; i < n; i++) { for (const auto& r: reaction_list()) { auto & rprop = m_graph[r].checked_property< Reaction >(); @@ -370,7 +411,9 @@ double Network::compute_all_reaction_rates(const unsigned n) const rprop.calc_rate(std::move(params)); } } - return get_time() - t_start; + clock_gettime(CLOCK_MONOTONIC_RAW, &t_end); + return TIME_DIFF(t_start, t_end)/5; + //return get_time() - t_start; } reaction_rate_t Network::get_reaction_rate(const Network::v_desc_t r) const @@ -766,6 +809,7 @@ void Network::print() const size_t num_inactive = 0ul; std::cout << "\n\nReactions:\n"; + double sum_rates = 0.0; for(const auto& vd : reaction_list()) { using directed_category = typename boost::graph_traits::directed_category; @@ -819,10 +863,14 @@ void Network::print() const std::cout << std::endl << " by the rate " << rp.get_rate() << " <= {" << rp.get_rate_formula() << "}" << std::endl; + sum_rates = sum_rates + rp.get_rate(); } // end of for loop over the reaction list std::cout << "Num inactive reactions: " << num_inactive << "/" << reaction_list().size() << std::endl; + std::cout << "Average of rates: " + << sum_rates << "/" << reaction_list().size() << ": " + << sum_rates / reaction_list().size() << std::endl; } /**@}*/ diff --git a/src/reaction_network/vertex.hpp b/src/reaction_network/vertex.hpp index d394bc4..06951dd 100644 --- a/src/reaction_network/vertex.hpp +++ b/src/reaction_network/vertex.hpp @@ -70,7 +70,7 @@ class Vertex { #if defined(WCS_HAS_SBML) template - Vertex(const LIBSBML_CPP_NAMESPACE::Species& species, const G& g); + Vertex(const LIBSBML_CPP_NAMESPACE::Model& model, const LIBSBML_CPP_NAMESPACE::Species& species, const G& g); template Vertex(const LIBSBML_CPP_NAMESPACE::Model& model, const @@ -153,7 +153,8 @@ Vertex::Vertex(const VertexFlat& flat, const G& g) #if defined(WCS_HAS_SBML) template -Vertex::Vertex(const LIBSBML_CPP_NAMESPACE::Species& species, const G& g) +Vertex::Vertex(const LIBSBML_CPP_NAMESPACE::Model& model, +const LIBSBML_CPP_NAMESPACE::Species& species, const G& g) : m_type(_species_), m_typeid(static_cast(_species_)), m_label(species.getIdAttribute()), @@ -162,15 +163,51 @@ Vertex::Vertex(const LIBSBML_CPP_NAMESPACE::Species& species, const G& g) { m_p = std::unique_ptr(new Species); - if (!isnan(species.getInitialAmount())) { + if (species.isSetInitialAmount()) { dynamic_cast(m_p.get())-> set_count(static_cast(species.getInitialAmount())); + // TODO: if units of species are mole then species Initial Amount has + // to be multiplied by Avogadro number + // TODO: species.getInitialConcentration() should be multiplied by // compartment volume and molarity (depending on the unit of // concentration) should be converted to count - } else if (!isnan(species.getInitialConcentration())) { + } else if (species.isSetInitialConcentration()) { + const LIBSBML_CPP_NAMESPACE::ListOfCompartments* compartment_list + = model.getListOfCompartments(); + const LIBSBML_CPP_NAMESPACE::ListOfUnitDefinitions* unit_definition_list + = model.getListOfUnitDefinitions(); + double compartment_size = compartment_list->get(species.getCompartment())->getSize(); + std::string compartment_unit ; + if (compartment_list->get(species.getCompartment())->isSetUnits()){ + compartment_unit = compartment_list->get(species.getCompartment())->getUnits(); + } else if (compartment_list->get(species.getCompartment())->getSpatialDimensions() == 3) { + compartment_unit = model.getVolumeUnits(); + } else if (compartment_list->get(species.getCompartment())->getSpatialDimensions() == 2) { + compartment_unit = model.getAreaUnits(); + } else if (compartment_list->get(species.getCompartment())->getSpatialDimensions() == 1) { + compartment_unit = model.getLengthUnits(); + } + const LIBSBML_CPP_NAMESPACE::UnitDefinition* unit_definition = unit_definition_list->get(compartment_unit); + const LIBSBML_CPP_NAMESPACE::ListOfUnits* unit_list + = unit_definition->getListOfUnits(); + unsigned int unitsSize = unit_list->size(); + double comp_unit = 1.0; + for (unsigned int iu = 0u; iu < unitsSize; iu++) { + const LIBSBML_CPP_NAMESPACE::Unit* unit = unit_list->get(iu); + comp_unit = comp_unit * pow(unit->getMultiplier()*pow(10,unit->getScale()),unit->getExponent()); + + } + double avog_num = 6.02214e+23; dynamic_cast(m_p.get())-> - set_count(static_cast(species.getInitialConcentration())); + // set_count(static_cast(species.getInitialConcentration() * (6.02214E23 * + // compartment_size) *comp_unit)); + set_count(static_cast(species.getInitialConcentration() * (6.02214E23 * + compartment_size) )); + // double avog_num = 6.02214e+23; + // printf("con %lf avo %lf\n", species.getInitialConcentration(), avog_num); + // while(1); + } } @@ -210,11 +247,31 @@ Vertex::Vertex( const LIBSBML_CPP_NAMESPACE::ListOfParameters* parameter_list = model.getListOfParameters(); unsigned int parametersSize = parameter_list->size(); - + //Add local parameters - const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list - = reaction.getKineticLaw()->getListOfLocalParameters(); - unsigned int num_local_parameters = local_parameter_list->size(); + if (model.getLevel() > 2) { + const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list + = reaction.getKineticLaw()->getListOfLocalParameters(); + unsigned int num_local_parameters = reaction.getKineticLaw()->getNumLocalParameters(); + + // Create an unordered_set for all local parameters of reaction + for (unsigned int pi = 0u; pi < num_local_parameters; pi++) { + const LIBSBML_CPP_NAMESPACE::LocalParameter* localparameter + = local_parameter_list->get(pi); + lpset.insert(localparameter->getIdAttribute()); + } + } else { + const LIBSBML_CPP_NAMESPACE::ListOfParameters* local_parameter_list + = reaction.getKineticLaw()->getListOfParameters(); + unsigned int num_local_parameters = reaction.getKineticLaw()->getNumParameters(); + // Create an unordered_set for all local parameters of reaction + for (unsigned int pi = 0u; pi < num_local_parameters; pi++) { + const LIBSBML_CPP_NAMESPACE::Parameter* localparameter + = local_parameter_list->get(pi); + lpset.insert(localparameter->getIdAttribute()); + } + } + sbml_utils sbml_o; pset=sbml_o.get_reaction_parameters(model, reaction); @@ -225,12 +282,6 @@ Vertex::Vertex( mpset.insert(parameter->getIdAttribute()); } - // Create an unordered_set for all local parameters of reaction - for (unsigned int pi = 0u; pi < num_local_parameters; pi++) { - const LIBSBML_CPP_NAMESPACE::LocalParameter* localparameter - = local_parameter_list->get(pi); - lpset.insert(localparameter->getIdAttribute()); - } // Put initial assignments in map using all_assignments_type @@ -284,16 +335,30 @@ Vertex::Vertex( } // reaction local parameters lpit = lpset.find(s_label) ; - if (lpit != lpset.end()) { - std::stringstream ss; - ss << local_parameter_list->get(s_label)->getValue(); - std::string parametervalue = ss.str(); - wholeformula = wholeformula + "var " + s_label - + " := " + parametervalue + "; "; - } + if (model.getLevel() > 2) { + const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list + = reaction.getKineticLaw()->getListOfLocalParameters(); + if (lpit != lpset.end()) { + std::stringstream ss; + ss << local_parameter_list->get(s_label)->getValue(); + std::string parametervalue = ss.str(); + wholeformula = wholeformula + "var " + s_label + + " := " + parametervalue + "; "; + } + } else { + const LIBSBML_CPP_NAMESPACE::ListOfParameters* local_parameter_list + = reaction.getKineticLaw()->getListOfParameters(); + if (lpit != lpset.end()) { + std::stringstream ss; + ss << local_parameter_list->get(s_label)->getValue(); + std::string parametervalue = ss.str(); + wholeformula = wholeformula + "var " + s_label + + " := " + parametervalue + "; "; + } + } } - //Add compartnents + //Add compartments const LIBSBML_CPP_NAMESPACE::ListOfCompartments* compartment_list = model.getListOfCompartments(); unsigned int compartmentsSize = compartment_list->size(); diff --git a/src/utils/generate_cxx_code.cpp b/src/utils/generate_cxx_code.cpp index 66f935e..ccb7222 100644 --- a/src/utils/generate_cxx_code.cpp +++ b/src/utils/generate_cxx_code.cpp @@ -170,6 +170,46 @@ generate_cxx_code::get_all_dependencies( return dependencies_no_dupl; } +static std::string add_braces_to_min( + char* math) +{ + std::string new_math = std::string(math); + std::string min_str ("min("); + std::string open_str ("("); + std::string close_str (")"); + std::size_t found = new_math.find(min_str); + while (found!=std::string::npos){ + new_math.insert(found, "std::"); + new_math.insert(found+9, 1 , '{'); + std::size_t close_pos = new_math.find(close_str,found+10); + std::size_t open_pos = new_math.find(open_str,found+10); + int num_open=1; + while( close_pos != std::string::npos) + { + if (open_pos == std::string::npos && num_open==1){ + new_math.insert(close_pos, 1 , '}'); + break; + } else if ( open_pos > close_pos ) { + if (num_open ==1) { + new_math.insert(close_pos, 1 , '}'); + break; + } else { + num_open = num_open - 1; + open_pos = new_math.find(open_str,close_pos+1); + close_pos = new_math.find(close_str,close_pos+1); + } + } else { + num_open = num_open + 1; + close_pos = new_math.find(close_str,open_pos+1); + open_pos = new_math.find(open_str,open_pos+1); + } + } + found = new_math.find(min_str,close_pos+1); + } + return new_math; +} + + void generate_cxx_code::return_denominators( const LIBSBML_CPP_NAMESPACE::ASTNode& formula, @@ -188,7 +228,8 @@ generate_cxx_code::return_denominators( } else { if ( std::find(denominators.cbegin(), denominators.cend(), SBML_formulaToString(denominator)) == denominators.cend()) { - denominators.push_back(SBML_formulaToString(denominator)); + std::string new_denominator = add_braces_to_min(SBML_formulaToString(denominator)); + denominators.push_back(new_denominator); } } } @@ -365,15 +406,29 @@ void generate_cxx_code::find_used_params( // Find used parameters in the rates for (unsigned int ic = 0u; ic < num_reactions; ic++) { const LIBSBML_CPP_NAMESPACE::Reaction& reaction = *(reaction_list->get(ic)); - const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list - = reaction.getKineticLaw()->getListOfLocalParameters(); - unsigned int num_localparameters = local_parameter_list->size(); - //genfile << "reaction" << reaction.getIdAttribute() <<": "; using reaction_parameters_t = std::unordered_set; reaction_parameters_t lpset; - for (unsigned int pi = 0u; pi < num_localparameters; pi++) { - const LIBSBML_CPP_NAMESPACE::LocalParameter* localparameter = local_parameter_list->get(pi); - lpset.insert(localparameter->getIdAttribute()); + // Check SBML level for local parameters + if (model.getLevel() > 2) { + const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list + = reaction.getKineticLaw()->getListOfLocalParameters(); + unsigned int num_localparameters = local_parameter_list->size(); + //genfile << "reaction" << reaction.getIdAttribute() <<": "; + + for (unsigned int pi = 0u; pi < num_localparameters; pi++) { + const LIBSBML_CPP_NAMESPACE::LocalParameter* localparameter = local_parameter_list->get(pi); + lpset.insert(localparameter->getIdAttribute()); + } + } else { + const LIBSBML_CPP_NAMESPACE::ListOfParameters* local_parameter_list + = reaction.getKineticLaw()->getListOfParameters(); + unsigned int num_localparameters = local_parameter_list->size(); + //genfile << "reaction" << reaction.getIdAttribute() <<": "; + + for (unsigned int pi = 0u; pi < num_localparameters; pi++) { + const LIBSBML_CPP_NAMESPACE::Parameter* localparameter = local_parameter_list->get(pi); + lpset.insert(localparameter->getIdAttribute()); + } } std::vector dependencies_set = get_all_dependencies(*reaction.getKineticLaw()->getMath(), @@ -1068,21 +1123,42 @@ void generate_cxx_code::print_reaction_rates( genfile << "#include \"" + header + '"' + "\n\n"; for (unsigned int ic = rid_start; ic < rid_end; ic++) { const LIBSBML_CPP_NAMESPACE::Reaction& reaction = *(reaction_list->get(ic)); - const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list - = reaction.getKineticLaw()->getListOfLocalParameters(); - unsigned int num_localparameters = local_parameter_list->size(); + unsigned int num_localparameters = 0u; + if (model.getLevel() > 2) { + const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list + = reaction.getKineticLaw()->getListOfLocalParameters(); + num_localparameters = local_parameter_list->size(); + } else { + const LIBSBML_CPP_NAMESPACE::ListOfParameters* local_parameter_list + = reaction.getKineticLaw()->getListOfParameters(); + num_localparameters = local_parameter_list->size(); + } genfile << "extern \"C\" " << Real << " wcs__rate_" << reaction.getIdAttribute() << "(const std::vector<" << Real <<">& __input) {\n"; - + + double sum_stoich = 0.0; //print reaction's local parameters using reaction_local_parameters_t = std::unordered_set; reaction_local_parameters_t localpset; - for (unsigned int pi = 0u; pi < num_localparameters; pi++) { - const LIBSBML_CPP_NAMESPACE::LocalParameter* localparameter = local_parameter_list->get(pi); - genfile << " constexpr " << Real << " " << localparameter->getIdAttribute() << " = " - << localparameter->getValue() << ";\n"; - localpset.insert(localparameter->getIdAttribute()); + if (model.getLevel() > 2) { + const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list + = reaction.getKineticLaw()->getListOfLocalParameters(); + for (unsigned int pi = 0u; pi < num_localparameters; pi++) { + const LIBSBML_CPP_NAMESPACE::LocalParameter* localparameter = local_parameter_list->get(pi); + //genfile << " constexpr " << Real << " " << localparameter->getIdAttribute() << " = " + // << localparameter->getValue() << ";\n"; + localpset.insert(localparameter->getIdAttribute()); + } + } else { + const LIBSBML_CPP_NAMESPACE::ListOfParameters* local_parameter_list + = reaction.getKineticLaw()->getListOfParameters(); + for (unsigned int pi = 0u; pi < num_localparameters; pi++) { + const LIBSBML_CPP_NAMESPACE::Parameter* localparameter = local_parameter_list->get(pi); + //genfile << " constexpr " << Real << " " << localparameter->getIdAttribute() << " = " + // << localparameter->getValue() << ";\n"; + localpset.insert(localparameter->getIdAttribute()); + } } //genfile << " int __ii=0;\n"; //genfile << "printf(\"Number of input arguments of %s received : %lu \", \"" @@ -1157,6 +1233,7 @@ void generate_cxx_code::print_reaction_rates( genfile << " wcs_global_var." << *it << " = wcs__rate_" << *it << "(" << function_input << ");\n"; } else { genfile << " " << Real << " " << *it << " = __input[" << par_index++ << "];\n"; + sum_stoich = sum_stoich + reaction.getReactant(*it)->getStoichiometry(); par_names.push_back(*it); } var_names.erase(*it); @@ -1199,7 +1276,172 @@ void generate_cxx_code::print_reaction_rates( par_names_nf.push_back(x); } } + const LIBSBML_CPP_NAMESPACE::ListOfUnitDefinitions* unit_definition_list + = model.getListOfUnitDefinitions(); + if (model.getLevel() > 2) { + int has_only_substance = 0; + double compartment_size = 0; + const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list + = reaction.getKineticLaw()->getListOfLocalParameters(); + const LIBSBML_CPP_NAMESPACE::ListOfCompartments* compartment_list + = model.getListOfCompartments(); + for (unsigned int pi = 0u; pi < num_localparameters; pi++) { + const LIBSBML_CPP_NAMESPACE::LocalParameter* localparameter = local_parameter_list->get(pi); + + // Check hasOnlySubstance + const ListOfSpecies* species_list = model.getListOfSpecies(); + const unsigned int num_species = species_list->size(); + for (unsigned int ic = 0u; ic < num_species; ic++) { + if (species_list->get(ic)->getHasOnlySubstanceUnits()) { + has_only_substance = 1; + } else { + compartment_size = compartment_list->get(species_list->get(ic)->getCompartment())->getSize(); + } + } + + if (localparameter->isSetUnits()){ + std::string local_par_unit = localparameter->getUnits(); + const LIBSBML_CPP_NAMESPACE::UnitDefinition* unit_definition = unit_definition_list->get(local_par_unit); + const LIBSBML_CPP_NAMESPACE::ListOfUnits* unit_list + = unit_definition->getListOfUnits(); + unsigned int unitsSize = unit_list->size(); + double comp_unit = 1.0; + for (unsigned int iu = 0u; iu < unitsSize; iu++) { + const LIBSBML_CPP_NAMESPACE::Unit* unit = unit_list->get(iu); + comp_unit = comp_unit * pow(unit->getMultiplier()*pow(10,unit->getScale()),unit->getExponent()); + } + genfile << " constexpr " << Real << " " << localparameter->getIdAttribute() << " = " + << localparameter->getValue() * comp_unit << ";\n"; + } else { + //Volume units/(substance units *timeunits) + if (model.isSetVolumeUnits() && model.isSetSubstanceUnits() && model.isSetTimeUnits()){ + std::string volume_units = model.getVolumeUnits(); + std::string substance_units = model.getSubstanceUnits(); + std::string time_units = model.getTimeUnits(); + double comp_unit = 1.0; + const LIBSBML_CPP_NAMESPACE::UnitDefinition* unit_definition; + const LIBSBML_CPP_NAMESPACE::ListOfUnits* unit_list; + unsigned int unitsSize = 0u; + + if (unit_definition_list->get(volume_units) != NULL) { + unit_definition = unit_definition_list->get(volume_units); + unit_list = unit_definition->getListOfUnits(); + unitsSize = unit_list->size(); + for (unsigned int iu = 0u; iu < unitsSize; iu++) { + const LIBSBML_CPP_NAMESPACE::Unit* unit = unit_list->get(iu); + comp_unit = comp_unit * pow(unit->getMultiplier()*pow(10,unit->getScale()),(sum_stoich-1) * unit->getExponent()); + } + } + + if (unit_definition_list->get(substance_units) != NULL) { + unit_definition = unit_definition_list->get(substance_units); + unit_list = unit_definition->getListOfUnits(); + unitsSize = unit_list->size(); + for (unsigned int iu = 0u; iu < unitsSize; iu++) { + const LIBSBML_CPP_NAMESPACE::Unit* unit = unit_list->get(iu); + if (unit->getScale() == -3) { + comp_unit = comp_unit * pow(unit->getMultiplier(),-((sum_stoich-1) * unit->getExponent())); + } else { + comp_unit = comp_unit * pow(unit->getMultiplier()*pow(10,unit->getScale()),-((sum_stoich-1) * unit->getExponent())); + } + } + } + + if (unit_definition_list->get(time_units) != NULL) { + unit_definition = unit_definition_list->get(time_units); + unit_list = unit_definition->getListOfUnits(); + unitsSize = unit_list->size(); + for (unsigned int iu = 0u; iu < unitsSize; iu++) { + const LIBSBML_CPP_NAMESPACE::Unit* unit = unit_list->get(iu); + comp_unit = comp_unit * pow(unit->getMultiplier()*pow(10,unit->getScale()),-((sum_stoich-1) * unit->getExponent())); + } + } + if (has_only_substance ==1) { + // genfile << " constexpr " << Real << " " << localparameter->getIdAttribute() << " = " + // << localparameter->getValue() * comp_unit << ";\n"; + genfile << " constexpr " << Real << " " << localparameter->getIdAttribute() << " = " + << localparameter->getValue() << ";\n"; + } else { + // genfile << " constexpr " << Real << " " << localparameter->getIdAttribute() << " = " + // << localparameter->getValue() * comp_unit / compartment_size << ";\n"; + // genfile << " constexpr " << Real << " " << localparameter->getIdAttribute() << " = " + // << localparameter->getValue() * comp_unit << ";\n"; + genfile << " constexpr " << Real << " " << localparameter->getIdAttribute() << " = " + << localparameter->getValue() << ";\n"; + } + } + } + } + } else { + const LIBSBML_CPP_NAMESPACE::ListOfParameters* local_parameter_list + = reaction.getKineticLaw()->getListOfParameters(); + for (unsigned int pi = 0u; pi < num_localparameters; pi++) { + const LIBSBML_CPP_NAMESPACE::Parameter* localparameter = local_parameter_list->get(pi); + + if (localparameter->isSetUnits()){ + std::string local_par_unit = localparameter->getUnits(); + const LIBSBML_CPP_NAMESPACE::UnitDefinition* unit_definition = unit_definition_list->get(local_par_unit); + const LIBSBML_CPP_NAMESPACE::ListOfUnits* unit_list + = unit_definition->getListOfUnits(); + unsigned int unitsSize = unit_list->size(); + double comp_unit = 1.0; + for (unsigned int iu = 0u; iu < unitsSize; iu++) { + const LIBSBML_CPP_NAMESPACE::Unit* unit = unit_list->get(iu); + comp_unit = comp_unit * pow(unit->getMultiplier()*pow(10,unit->getScale()),unit->getExponent()); + } + genfile << " constexpr " << Real << " " << localparameter->getIdAttribute() << " = " + << localparameter->getValue() * comp_unit << ";\n"; + } else { + //Volume units/(substance units *timeunits) + if (model.isSetVolumeUnits() && model.isSetSubstanceUnits() && model.isSetTimeUnits()){ + std::string volume_units = model.getVolumeUnits(); + std::string substance_units = model.getSubstanceUnits(); + std::string time_units = model.getTimeUnits(); + double comp_unit = 1.0; + const LIBSBML_CPP_NAMESPACE::UnitDefinition* unit_definition; + const LIBSBML_CPP_NAMESPACE::ListOfUnits* unit_list; + unsigned int unitsSize = 0u; + + if (unit_definition_list->get(volume_units) != NULL) { + unit_definition = unit_definition_list->get(volume_units); + unit_list = unit_definition->getListOfUnits(); + unitsSize = unit_list->size(); + for (unsigned int iu = 0u; iu < unitsSize; iu++) { + const LIBSBML_CPP_NAMESPACE::Unit* unit = unit_list->get(iu); + comp_unit = comp_unit * pow(unit->getMultiplier()*pow(10,unit->getScale()),(sum_stoich-1) * unit->getExponent()); + } + } + + if (unit_definition_list->get(substance_units) != NULL) { + unit_definition = unit_definition_list->get(substance_units); + unit_list = unit_definition->getListOfUnits(); + unitsSize = unit_list->size(); + for (unsigned int iu = 0u; iu < unitsSize; iu++) { + const LIBSBML_CPP_NAMESPACE::Unit* unit = unit_list->get(iu); + if (unit->getScale() == -3) { + comp_unit = comp_unit * pow(unit->getMultiplier(),-((sum_stoich-1) * unit->getExponent())); + } else { + comp_unit = comp_unit * pow(unit->getMultiplier()*pow(10,unit->getScale()),-((sum_stoich-1) * unit->getExponent())); + } + } + } + if (unit_definition_list->get(time_units) != NULL) { + unit_definition = unit_definition_list->get(time_units); + unit_list = unit_definition->getListOfUnits(); + unitsSize = unit_list->size(); + for (unsigned int iu = 0u; iu < unitsSize; iu++) { + const LIBSBML_CPP_NAMESPACE::Unit* unit = unit_list->get(iu); + comp_unit = comp_unit * pow(unit->getMultiplier()*pow(10,unit->getScale()),-((sum_stoich-1) * unit->getExponent())); + } + } + + genfile << " constexpr " << Real << " " << localparameter->getIdAttribute() << " = " + << localparameter->getValue() * comp_unit << ";\n"; + } + } + } + } // put input parameters into maps dep_params_f.insert(std::make_pair(reaction.getIdAttribute(), par_names)); @@ -1232,8 +1474,9 @@ void generate_cxx_code::print_reaction_rates( } math = *reaction.getKineticLaw()->getMath(); update_scope_ast_node(math, wcs_all_var, wcs_all_const, localpset); + std::string new_math = add_braces_to_min(SBML_formulaToString(&math)); genfile << " " << Real << " " << reaction.getIdAttribute() << " = " - << SBML_formulaToString( &math) + << new_math //SBML_formulaToString( &math) << ";\n"; genfile << " if (!isfinite(" << reaction.getIdAttribute() << ")) {\n"; // find denominators of the dependent expressions of the reaction rate formula @@ -1674,6 +1917,7 @@ void generate_cxx_code::write_header( << "#include \n" << "#include \n" << "#include \n" + << "#include \n" << "#include \"utils/exception.hpp\"\n\n" << "//Include the text of the SBML in here.\n" << "//Use \"xxd -i\" to get the text as a C symbol.\n" diff --git a/src/utils/graph_factory.hpp b/src/utils/graph_factory.hpp index 46e8b4f..b6a7504 100644 --- a/src/utils/graph_factory.hpp +++ b/src/utils/graph_factory.hpp @@ -303,7 +303,6 @@ GraphFactory::convert_to( { g.m_vertices[vd].m_in_edges.reserve(num_reactants); } - // Add reactants species for (unsigned int si = 0u; si < num_reactants; si++) { const auto &reactant = *(reaction.getReactant(si)); @@ -315,7 +314,8 @@ GraphFactory::convert_to( reaction_in_species.insert(s_label); if (it == smap.cend()) { - wcs::Vertex vs(*species_list->get(reactant.getSpecies()), g); + wcs::Vertex vs(model, *species_list->get(reactant.getSpecies()), g); + //WCS_THROW("TESTT"); vds = boost::add_vertex(vs, g); smap.insert(std::make_pair(s_label,vds)); } else { @@ -323,27 +323,38 @@ GraphFactory::convert_to( } std::string e_label = g[vds].get_label() + '|' + g[vd].get_label(); - eit = emap.find(e_label); + if (!species_list->get(reactant.getSpecies())->getConstant()) { //Add as a reactant + eit = emap.find(e_label); + + if (eit == emap.cend()) { + const auto ret = boost::add_edge(vds, vd, g); - if (eit == emap.cend()) { + if (!ret.second) { + WCS_THROW("Please check the reactions in your SBML file"); + return; + } + g[ret.first].set_stoichiometry_ratio(reactant.getStoichiometry()); + g[ret.first].set_label(e_label); + emap.insert(std::make_pair(e_label, ret.first)); + } else { + const auto& edge_found = eit->second; + stoic_t new_stoich = g[edge_found].get_stoichiometry_ratio() + + reactant.getStoichiometry(); + g[edge_found].set_stoichiometry_ratio(new_stoich); + } + } else { //Add as a modifier const auto ret = boost::add_edge(vds, vd, g); if (!ret.second) { WCS_THROW("Please check the reactions in your SBML file"); return; } - g[ret.first].set_stoichiometry_ratio(reactant.getStoichiometry()); + g[ret.first].set_stoichiometry_ratio(0); g[ret.first].set_label(e_label); emap.insert(std::make_pair(e_label, ret.first)); - } else { - const auto& edge_found = eit->second; - stoic_t new_stoich = g[edge_found].get_stoichiometry_ratio() - + reactant.getStoichiometry(); - g[edge_found].set_stoichiometry_ratio(new_stoich); - } + } } - // Add modifiers species for (unsigned int si = 0u; si < num_modifiers; si++) { const auto &modifier = *(reaction.getModifier(si)); @@ -355,7 +366,7 @@ GraphFactory::convert_to( reaction_in_species.insert(s_label); if (it == smap.cend()) { - wcs::Vertex vs(*species_list->get(modifier.getSpecies()), g); + wcs::Vertex vs(model, *species_list->get(modifier.getSpecies()), g); vds = boost::add_vertex(vs, g); smap.insert(std::make_pair(s_label,vds)); } else { @@ -396,7 +407,7 @@ GraphFactory::convert_to( reaction_in_species.insert(s_label); if (it == smap.cend()) { - wcs::Vertex vs(*species_list->get(x), g); + wcs::Vertex vs(model, *species_list->get(x), g); vds = boost::add_vertex(vs, g); smap.insert(std::make_pair(s_label,vds)); } else { @@ -434,7 +445,7 @@ GraphFactory::convert_to( if (risit == reaction_in_species.cend()) { reaction_in_species.insert(x); if (it == smap.cend()) { - wcs::Vertex vs(*species_list->get(x), g); + wcs::Vertex vs(model, *species_list->get(x), g); vds = boost::add_vertex(vs, g); smap.insert(std::make_pair(x,vds)); } else { @@ -472,7 +483,7 @@ GraphFactory::convert_to( risit = reaction_in_species.find(x); if (risit == reaction_in_species.cend()) { if (it == smap.cend()) { - wcs::Vertex vs(*species_list->get(x), g); + wcs::Vertex vs(model, *species_list->get(x), g); vds = boost::add_vertex(vs, g); smap.insert(std::make_pair(x,vds)); } else { @@ -506,7 +517,7 @@ GraphFactory::convert_to( v_new_desc_t vds; if (it == smap.cend()) { - wcs::Vertex vs(*species_list->get(product.getSpecies()), g); + wcs::Vertex vs(model, *species_list->get(product.getSpecies()), g); vds = boost::add_vertex(vs, g); smap.insert(std::make_pair(s_label,vds)); } else { @@ -514,22 +525,35 @@ GraphFactory::convert_to( } std::string e_label = g[vd].get_label() + '|' + g[vds].get_label(); - eit = emap.find(e_label); - if (eit == emap.cend()) { - const auto ret = boost::add_edge(vd, vds, g); + + if (!species_list->get(product.getSpecies())->getConstant()) { //Add as a product + eit = emap.find(e_label); + if (eit == emap.cend()) { + const auto ret = boost::add_edge(vd, vds, g); + + if (!ret.second) { + WCS_THROW("Please check the reactions in your SBML file"); + return; + } + g[ret.first].set_stoichiometry_ratio(product.getStoichiometry()); + g[ret.first].set_label( g[vd].get_label() + '|' + g[vds].get_label()); + emap.insert(std::make_pair(e_label, ret.first)); + } else { + const auto& edge_found = eit->second; + stoic_t new_stoich = g[edge_found].get_stoichiometry_ratio() + + product.getStoichiometry(); + g[edge_found].set_stoichiometry_ratio(new_stoich); + } + } else { //Add as a modifier + const auto ret = boost::add_edge(vds, vd, g); if (!ret.second) { WCS_THROW("Please check the reactions in your SBML file"); return; } - g[ret.first].set_stoichiometry_ratio(product.getStoichiometry()); - g[ret.first].set_label( g[vd].get_label() + '|' + g[vds].get_label()); + g[ret.first].set_stoichiometry_ratio(0); + g[ret.first].set_label(e_label); emap.insert(std::make_pair(e_label, ret.first)); - } else { - const auto& edge_found = eit->second; - stoic_t new_stoich = g[edge_found].get_stoichiometry_ratio() - + product.getStoichiometry(); - g[edge_found].set_stoichiometry_ratio(new_stoich); } } } diff --git a/src/utils/sbml_utils.cpp b/src/utils/sbml_utils.cpp index 9efdba7..7599cca 100644 --- a/src/utils/sbml_utils.cpp +++ b/src/utils/sbml_utils.cpp @@ -81,12 +81,22 @@ sbml_utils::find_undeclared_species_in_reaction_formula( } // create an unordered_set for reaction local parameters - const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list - = reaction.getKineticLaw()->getListOfLocalParameters(); - unsigned int num_local_parameters = local_parameter_list->size(); + if (model.getLevel() > 2) { + const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list + = reaction.getKineticLaw()->getListOfLocalParameters(); + unsigned int num_local_parameters = local_parameter_list->size(); - for (unsigned int pi = 0u; pi < num_local_parameters; pi++) { - lpset.insert(local_parameter_list->get(pi)->getIdAttribute()); + for (unsigned int pi = 0u; pi < num_local_parameters; pi++) { + lpset.insert(local_parameter_list->get(pi)->getIdAttribute()); + } + } else { + const LIBSBML_CPP_NAMESPACE::ListOfParameters* local_parameter_list + = reaction.getKineticLaw()->getListOfParameters(); + unsigned int num_local_parameters = local_parameter_list->size(); + + for (unsigned int pi = 0u; pi < num_local_parameters; pi++) { + lpset.insert(local_parameter_list->get(pi)->getIdAttribute()); + } } @@ -94,7 +104,11 @@ sbml_utils::find_undeclared_species_in_reaction_formula( unsigned int num_reactants = reaction.getNumReactants(); for (unsigned int ri = 0u; ri < num_reactants; ri++) { const auto &reactant = *(reaction.getReactant(ri)); - rset.insert(species_list->get(reactant.getSpecies())->getIdAttribute()); + if (!species_list->get(reactant.getSpecies())->getConstant()) { //Add as a reactant + rset.insert(species_list->get(reactant.getSpecies())->getIdAttribute()); + } else { + mset.insert(species_list->get(reactant.getSpecies())->getIdAttribute()); + } } // create an unordered set for reaction_modifiers @@ -174,7 +188,9 @@ sbml_utils::get_reaction_parameters( for (unsigned int ri = 0u; ri < num_reactants; ri++) { const auto &reactant = *(reaction.getReactant(ri)); - rset.insert(species_list->get(reactant.getSpecies())->getIdAttribute()); + if (!species_list->get(reactant.getSpecies())->getConstant()) { //Add as a reactant + rset.insert(species_list->get(reactant.getSpecies())->getIdAttribute()); + } } // create an unordered_set for model compartments From a2edf786fec0a4d9731199d4eb190461b6d24441 Mon Sep 17 00:00:00 2001 From: Konstantia Georgouli Date: Mon, 7 Jun 2021 10:33:22 -0700 Subject: [PATCH 02/14] Fix min mathematical function and support models with initial concentration --- src/reaction_network/network.cpp | 66 +---- src/reaction_network/vertex.cpp | 15 +- src/reaction_network/vertex.hpp | 21 +- src/utils/CMakeLists.txt | 2 + src/utils/FunctionInterfaceForJIT.cpp | 343 ++++++++++++++++++++++++++ src/utils/FunctionInterfaceForJIT.hpp | 74 ++++++ src/utils/generate_cxx_code.cpp | 332 +++++++++---------------- 7 files changed, 575 insertions(+), 278 deletions(-) create mode 100644 src/utils/FunctionInterfaceForJIT.cpp create mode 100644 src/utils/FunctionInterfaceForJIT.hpp diff --git a/src/reaction_network/network.cpp b/src/reaction_network/network.cpp index 43eca39..b6d504a 100644 --- a/src/reaction_network/network.cpp +++ b/src/reaction_network/network.cpp @@ -22,11 +22,6 @@ #include // is_same<> #include // lexicographical_compare(), sort() #include // numeric_limits -#include - -#define TIME_DIFF(Tstart, Tend ) \ - ((double) (1000000000L * ((Tend).tv_sec - (Tstart).tv_sec) + \ - (Tend).tv_nsec - (Tstart).tv_nsec) / 1000000000L) #if defined(WCS_HAS_SBML) #include @@ -160,12 +155,8 @@ void Network::loadSBML(const std::string sbml_filename, const bool reuse) typename params_map_t::const_iterator pit; struct timespec t_1, t_2; - clock_gettime(CLOCK_MONOTONIC_RAW, &t_1); code_generator.generate_code(*model, m_dep_params_f, m_dep_params_nf, m_rate_rules_dep_map); - clock_gettime(CLOCK_MONOTONIC_RAW, &t_2); - using std::operator<<; - std::cout << "Time for creating the generated file: " << TIME_DIFF(t_1, t_2 ) << " sec" << std::endl; const std::string library_file = code_generator.compile_code(); using std::operator<<; @@ -201,8 +192,6 @@ void Network::init() m_species.reserve(num_vertices); v_iter_t vi, vi_end; - int reaction_sz = 0; - double sum_in = 0, sum_reac_mod = 0, sum_involved = 0, sum_out=0; for (boost::tie(vi, vi_end) = boost::vertices(m_graph); vi != vi_end; ++vi) { const v_prop_t& v = m_graph[*vi]; @@ -210,8 +199,6 @@ void Network::init() if (vt == v_prop_t::_species_) { m_species.emplace_back(*vi); } else { - reaction_sz = reaction_sz +1; - int reaction_in = 0, reac_mod_in = 0, reaction_out =0; using directed_category = boost::graph_traits::directed_category; constexpr bool is_bidirectional = std::is_same::value; @@ -241,22 +228,18 @@ void Network::init() boost::make_iterator_range(boost::in_edges(reaction, m_graph))) { v_desc_t reactant = boost::source(ei_in, m_graph); - // #if !defined(WCS_HAS_EXPRTK) - // //check for the reactants which are not actually reactants and put their stoichiometry 0 - // if ( std::find(params_reactants.begin(), params_reactants.end(), - // m_graph[reactant].get_label()) == params_reactants.end()) { - // //std::cout << "Not reactant " << m_graph[reactant].get_label() << std::endl; - // m_graph[ei_in].set_stoichiometry_ratio(0); - // } - // #endif // !defined(WCS_HAS_EXPRTK) + #if !defined(WCS_HAS_EXPRTK) + //check for the reactants which are not actually reactants and put their stoichiometry 0 + if ( std::find(params_reactants.begin(), params_reactants.end(), + m_graph[reactant].get_label()) == params_reactants.end()) { + //std::cout << "Not reactant " << m_graph[reactant].get_label() << std::endl; + m_graph[ei_in].set_stoichiometry_ratio(0); + } + #endif // !defined(WCS_HAS_EXPRTK) const auto st = m_graph[ei_in].get_stoichiometry_ratio(); involved_species.insert(std::make_pair(m_graph[reactant].get_label(), std::make_pair(reactant, st))); - if (m_graph[ei_in].get_stoichiometry_ratio() > 0 ){ - reaction_in = reaction_in + 1; - } - reac_mod_in = reac_mod_in + 1; } } @@ -265,7 +248,6 @@ void Network::init() v_desc_t product = boost::target(ei_out, m_graph); products.insert(std::make_pair(m_graph[product].get_label(), std::make_pair(product, 1))); - reaction_out = reaction_out + 1; } m_reactions.emplace_back(reaction); @@ -315,31 +297,12 @@ void Network::init() r.set_rate_inputs(involved_species); #endif // !defined(WCS_HAS_EXPRTK) - auto & rprop = m_graph[*vi].checked_property< Reaction >(); - const auto& ri = rprop.get_rate_inputs(); - - sum_involved = sum_involved + ri.size(); - sum_in = sum_in + reaction_in; - sum_reac_mod = sum_reac_mod + reac_mod_in; - sum_out = sum_out + reaction_out; - //using std::operator<<; - //std::cout << "involved species: " << ri.size() << " reactants:" << reaction_in << std::endl; - r.set_products(products); set_reaction_rate(*vi); } } - double mean_in = sum_in / reaction_sz; - double mean_out = sum_out / reaction_sz; - double mean_reac_mod = sum_reac_mod / reaction_sz; - double mean_involved = sum_involved / reaction_sz; - using std::operator<<; - - std::cout << "Reactions " << reaction_sz << " mean of inputs " << mean_reac_mod - << ", mean of reactants " << mean_in << ", mean of species involved in reaction rate " << mean_involved - << " and, mean of products " << mean_out << std::endl; sort_species(); build_index_maps(); @@ -390,9 +353,7 @@ reaction_rate_t Network::set_reaction_rate(const Network::v_desc_t r) const double Network::compute_all_reaction_rates(const unsigned n) const { - //double t_start = get_time(); - struct timespec t_start, t_end; - clock_gettime(CLOCK_MONOTONIC_RAW, &t_start); + double t_start = get_time(); for (unsigned i = 0u; i < n; i++) { for (const auto& r: reaction_list()) { auto & rprop = m_graph[r].checked_property< Reaction >(); @@ -411,9 +372,7 @@ double Network::compute_all_reaction_rates(const unsigned n) const rprop.calc_rate(std::move(params)); } } - clock_gettime(CLOCK_MONOTONIC_RAW, &t_end); - return TIME_DIFF(t_start, t_end)/5; - //return get_time() - t_start; + return get_time() - t_start; } reaction_rate_t Network::get_reaction_rate(const Network::v_desc_t r) const @@ -809,7 +768,6 @@ void Network::print() const size_t num_inactive = 0ul; std::cout << "\n\nReactions:\n"; - double sum_rates = 0.0; for(const auto& vd : reaction_list()) { using directed_category = typename boost::graph_traits::directed_category; @@ -863,14 +821,10 @@ void Network::print() const std::cout << std::endl << " by the rate " << rp.get_rate() << " <= {" << rp.get_rate_formula() << "}" << std::endl; - sum_rates = sum_rates + rp.get_rate(); } // end of for loop over the reaction list std::cout << "Num inactive reactions: " << num_inactive << "/" << reaction_list().size() << std::endl; - std::cout << "Average of rates: " - << sum_rates << "/" << reaction_list().size() << ": " - << sum_rates / reaction_list().size() << std::endl; } /**@}*/ diff --git a/src/reaction_network/vertex.cpp b/src/reaction_network/vertex.cpp index 859c83f..df62d22 100644 --- a/src/reaction_network/vertex.cpp +++ b/src/reaction_network/vertex.cpp @@ -30,7 +30,7 @@ std::map Vertex::str_vt { }; Vertex::Vertex() -: m_type(_undefined_), m_label(""), +: m_type(_undefined_), m_label(""), m_annotation(""), m_pid(unassigned_partition), m_p(nullptr) { @@ -41,6 +41,7 @@ Vertex::Vertex(const Vertex& rhs) : m_type(rhs.m_type), m_typeid(rhs.m_typeid), m_label(rhs.m_label), + m_annotation(rhs.m_annotation), m_pid(rhs.m_pid), m_p((!rhs.m_p)? nullptr : rhs.m_p.get()->clone()) { @@ -53,6 +54,7 @@ Vertex::Vertex(Vertex&& rhs) noexcept { if (this != &rhs) { m_label = std::move(rhs.m_label); + m_annotation = std::move(rhs.m_annotation); m_p = std::move(rhs.m_p); reset(rhs); @@ -65,6 +67,7 @@ Vertex& Vertex::operator=(const Vertex& rhs) m_type = rhs.m_type; m_typeid = rhs.m_typeid; m_label = rhs.m_label; + m_annotation = rhs.m_annotation; m_pid = rhs.m_pid; m_p = (!rhs.m_p)? nullptr : rhs.m_p.get()->clone(); } @@ -77,6 +80,7 @@ Vertex& Vertex::operator=(Vertex&& rhs) noexcept m_type = rhs.m_type; m_typeid = rhs.m_typeid; m_label = std::move(rhs.m_label); + m_annotation = std::move(rhs.m_annotation); m_pid = rhs.m_pid; m_p = std::move(rhs.m_p); @@ -150,6 +154,15 @@ partition_id_t Vertex::get_partition() const { return m_pid; } +void Vertex::set_annotation(const std::string& f) +{ + m_annotation = f; +} + +std::string Vertex::get_annotation() const +{ + return m_annotation; +} /**@}*/ } // end of namespace wcs diff --git a/src/reaction_network/vertex.hpp b/src/reaction_network/vertex.hpp index 06951dd..4549cd6 100644 --- a/src/reaction_network/vertex.hpp +++ b/src/reaction_network/vertex.hpp @@ -89,6 +89,8 @@ class Vertex { std::string get_label() const; void set_partition(const partition_id_t pid); partition_id_t get_partition() const; + void set_annotation(const std::string& f); + std::string get_annotation() const; template P& property() const; template P& checked_property() const; @@ -106,6 +108,7 @@ class Vertex { int m_typeid; ///< The vertex type in integer form used for GraphML parsing std::string m_label; ///< The vertex label partition_id_t m_pid; ///< The id of the partition to which this edge belongs + std::string m_annotation; ///< vertex annotation /** * The pointer to the detailed property object, which is polymorphic. @@ -196,18 +199,11 @@ const LIBSBML_CPP_NAMESPACE::Species& species, const G& g) for (unsigned int iu = 0u; iu < unitsSize; iu++) { const LIBSBML_CPP_NAMESPACE::Unit* unit = unit_list->get(iu); comp_unit = comp_unit * pow(unit->getMultiplier()*pow(10,unit->getScale()),unit->getExponent()); - } double avog_num = 6.02214e+23; dynamic_cast(m_p.get())-> - // set_count(static_cast(species.getInitialConcentration() * (6.02214E23 * - // compartment_size) *comp_unit)); set_count(static_cast(species.getInitialConcentration() * (6.02214E23 * - compartment_size) )); - // double avog_num = 6.02214e+23; - // printf("con %lf avo %lf\n", species.getInitialConcentration(), avog_num); - // while(1); - + compartment_size) *comp_unit)); } } @@ -246,8 +242,7 @@ Vertex::Vertex( //Add parameters const LIBSBML_CPP_NAMESPACE::ListOfParameters* parameter_list = model.getListOfParameters(); - unsigned int parametersSize = parameter_list->size(); - + unsigned int parametersSize = parameter_list->size(); //Add local parameters if (model.getLevel() > 2) { const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list @@ -382,6 +377,12 @@ Vertex::Vertex( wholeformula = wholeformula + "m_rate := " + formula + ";"; dynamic_cast*>(m_p.get())->set_rate_formula(wholeformula); + if (reaction.isSetAnnotation()) { + std::string annot= reaction.getAnnotationString() ; + m_annotation = annot; + } else { + m_annotation = ""; + } #if !defined(WCS_HAS_EXPRTK) dynamic_cast*>(m_p.get())->ReactionBase::set_calc_rate_fn(reaction_function_rate); diff --git a/src/utils/CMakeLists.txt b/src/utils/CMakeLists.txt index 47ccdb9..2faac5a 100644 --- a/src/utils/CMakeLists.txt +++ b/src/utils/CMakeLists.txt @@ -12,6 +12,7 @@ set_full_path(THIS_DIR_HEADERS rngen_impl.hpp samples_ssa.hpp sbml_utils.hpp + FunctionInterfaceForJIT.hpp seed.hpp state_io.hpp state_io_impl.hpp @@ -39,6 +40,7 @@ set_full_path(THIS_DIR_SOURCES print_vertices.cpp samples_ssa.cpp sbml_utils.cpp + FunctionInterfaceForJIT.cpp trajectory.cpp trace_ssa.cpp trace_generic.cpp diff --git a/src/utils/FunctionInterfaceForJIT.cpp b/src/utils/FunctionInterfaceForJIT.cpp new file mode 100644 index 0000000..b85f2be --- /dev/null +++ b/src/utils/FunctionInterfaceForJIT.cpp @@ -0,0 +1,343 @@ +/****************************************************************************** + * * + * 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 * + * * + ******************************************************************************/ + +#include "utils/FunctionInterfaceForJIT.hpp" + +#if defined(WCS_HAS_CONFIG) +#include "wcs_config.hpp" +#else +#error "no config" +#endif + + +#if defined(WCS_HAS_SBML) +#include +#include +#endif // defined(WCS_HAS_SBML) + +#include +#include +// #include +// #include + +namespace wcs { +/** \addtogroup wcs_utils + * @{ */ + +FunctionInterfaceForJIT::FunctionInterfaceForJIT() +{} + +#if defined(WCS_HAS_SBML) +using ast_function_ptr = std::function; +using node_structure = std::tuple; +using node_map = std::unordered_map; +node_map nodetypes = { + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_ABS,{nullptr,"std::abs(", "", ")"}}, + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_ARCCOS,{nullptr,"std::acos(", "", ")"}}, + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_ARCCOSH,{nullptr,"std::acosh(", "", ")"}}, + // {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_ARCCOT,{nullptr,"std::(", "", ")"}}, // arccotangent + // {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_ARCCOTH,{nullptr,"std::(", "", ")"}}, // Hyperbolic arccotangent + // {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_ARCCSC,{nullptr,"std::(", "", ")"}}, // arccosecant + // {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_ARCCSCH,{nullptr,"std::(", "", ")"}}, // Hyperbolic arccosecant + // {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_ARCSEC,{nullptr,"std::(", "", ")"}}, // arcsecant + // {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_ARCSECH,{nullptr,"std::(", "", ")"}}, // Hyperbolic arcsecant + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_ARCSIN,{nullptr,"std::asin(", "", ")"}}, + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_ARCSINH,{nullptr,"std::asinh(", "", ")"}}, + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_ARCTAN,{nullptr,"std::atan(", "", ")"}}, + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_ARCTANH,{nullptr,"std::atanh(", "", ")"}}, + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_CEILING,{nullptr,"std::ceil(", "", ")"}}, + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_COS,{nullptr,"std::cos(", "", ")"}}, + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_COSH,{nullptr,"std::cosh(", "", ")"}}, + // {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_COT,{nullptr,"std::(", "", ")"}}, //cotangent + // {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_COTH,{nullptr,"std::(", "", ")"}}, // Hyperbolic cotangent + // {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_CSC,{nullptr,"std::(", "", ")"}}, // cosecant + // {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_CSCH,{nullptr,"std::(", "", ")"}}, // Hyperbolic cosecant + // {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_DELAY,{nullptr,"std::(", "", ")"}}, //Delay + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_EXP,{nullptr,"std::exp(", "", ")"}}, + // {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_FACTORIAL,{nullptr,"std::(", "", ")"}}, // Factorial + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_FLOOR,{nullptr,"std::floor(", "", ")"}}, + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_LN,{nullptr,"std::log(", "", ")"}}, // natural log + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_LOG,{nullptr,"", "", ""}}, // function pointer------- + // {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_PIECEWISE,{nullptr,"std::(", "", ")"}}, // Piecewise + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_POWER,{nullptr,"std::pow(", ", ", ")"}}, // NOT SURE + // {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_SEC,{nullptr,"std::(", "", ")"}}, // Secant + // {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_SECH,{nullptr,"std::(", "", ")"}}, // Hyperbolic Secant + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_SIN,{nullptr,"std::sin(", "", ")"}}, + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_SINH,{nullptr,"std::sinh(", "", ")"}}, + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_TAN,{nullptr,"std::tan(", "", ")"}}, + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_TANH,{nullptr,"std::tanh(", "", ")"}}, + + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_MIN,{nullptr,"std::min({", ", ", "})"}}, + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_MAX,{nullptr,"std::max({", ", ", "})"}}, + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_QUOTIENT,{nullptr,"", "", ""}}, // function pointer div ------- + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_RATE_OF,{nullptr,"", "", ""}}, // function pointer rate_of ------- + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_REM,{nullptr,"std::remainder(", ", ", ")"}}, + + {LIBSBML_CPP_NAMESPACE::AST_FUNCTION_ROOT,{nullptr,"", "", ""}}, //function pointer + }; +/** + * Visits the given ASTNode as a function. For this node only the + * traversal is preorder. + */ +/* WCS addition: Add "std::" before min and "{...}" between parenthesis in order to support min */ +void FunctionInterfaceForJIT::FormulaFormatter_visitFunction_wcs ( const ASTNode_t *parent, + const ASTNode_t *node, + strbuff_t *sb ) +{ + unsigned int numChildren = ASTNode_getNumChildren(node); + unsigned int n; + + // int node_type = ASTNode_getType(node); + LIBSBML_CPP_NAMESPACE::ASTNodeType_t node_type = ASTNode_getType(node); + constexpr LIBSBML_CPP_NAMESPACE::ASTNodeType_t ast_func_min = LIBSBML_CPP_NAMESPACE::AST_FUNCTION_MIN; + constexpr LIBSBML_CPP_NAMESPACE::ASTNodeType_t ast_func_max = LIBSBML_CPP_NAMESPACE::AST_FUNCTION_MAX; + if (node_type == ast_func_max || node_type == ast_func_min) { + StringBuffer_append(sb, "std::"); + } + FormulaFormatter_format(sb, node); + + // support min in sbml + StringBuffer_appendChar(sb, '('); + if (node_type == 321 || node_type == 320) { + StringBuffer_appendChar(sb, '{'); + } + + if (numChildren > 0) + { + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs( node, LIBSBML_CPP_NAMESPACE::ASTNode_getChild(node, 0), sb ); + } + + for (n = 1; n < numChildren; n++) + { + StringBuffer_appendChar(sb, ','); + StringBuffer_appendChar(sb, ' '); + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs( node, LIBSBML_CPP_NAMESPACE::ASTNode_getChild(node, n), sb ); + } + + if (node_type == 321 || node_type == 320) { + StringBuffer_appendChar(sb, '}'); + } + StringBuffer_appendChar(sb, ')'); +} + + +/** + * Visits the given ASTNode as the function "log(10, x)" and in doing so, + * formats it as "log10(x)" (where x is any subexpression). + */ +/* WCS addition: Change function call from FormulaFormatter_visit to FormulaFormatter_visit_wcs */ +void FunctionInterfaceForJIT::FormulaFormatter_visitLog10_wcs ( const ASTNode_t *parent, + const ASTNode_t *node, + strbuff_t *sb ) +{ + StringBuffer_append(sb, "log10("); + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs(node, ASTNode_getChild(node, 1), sb); + StringBuffer_appendChar(sb, ')'); +} + + +/** + * Visits the given ASTNode as the function "root(2, x)" and in doing so, + * formats it as "sqrt(x)" (where x is any subexpression). + */ +/* WCS addition: Change function call from FormulaFormatter_visit to FormulaFormatter_visit_wcs */ +void FunctionInterfaceForJIT::FormulaFormatter_visitSqrt_wcs ( const ASTNode_t *parent, + const ASTNode_t *node, + strbuff_t *sb ) +{ + StringBuffer_append(sb, "sqrt("); + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs(node, ASTNode_getChild(node, 1), sb); + StringBuffer_appendChar(sb, ')'); +} + + +/** + * Visits the given ASTNode as a unary minus. For this node only the + * traversal is preorder. + */ +/* WCS addition: Change function call from FormulaFormatter_visit to FormulaFormatter_visit_wcs */ +void FunctionInterfaceForJIT::FormulaFormatter_visitUMinus_wcs ( const ASTNode_t *parent, + const ASTNode_t *node, + strbuff_t *sb ) +{ + StringBuffer_appendChar(sb, '-'); + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs ( node, ASTNode_getLeftChild(node), sb ); +} + + +/** + * Visits the given ASTNode and continues the inorder traversal. + */ +/* WCS addition: Change function call from FormulaFormatter_visit to FormulaFormatter_visit_wcs */ +void FunctionInterfaceForJIT::FormulaFormatter_visitOther_wcs ( const ASTNode_t *parent, + const ASTNode_t *node, + strbuff_t *sb ) +{ + unsigned int numChildren = ASTNode_getNumChildren(node); + int group = FormulaFormatter_isGrouped(parent, node); + unsigned int n; + + + if (group) + { + StringBuffer_appendChar(sb, '('); + } + + if (numChildren == 0) { + FormulaFormatter_format(sb, node); + } + + else if (numChildren == 1) + { + //I believe this would only be called for invalid ASTNode setups, + // but this could in theory occur. This is the safest + // behavior I can think of. + FormulaFormatter_format(sb, node); + StringBuffer_appendChar(sb, '('); + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs( node, ASTNode_getChild(node, 0), sb ); + StringBuffer_appendChar(sb, ')'); + } + + else { + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs( node, ASTNode_getChild(node, 0), sb ); + + for (n = 1; n < numChildren; n++) + { + FormulaFormatter_format(sb, node); + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs( node, ASTNode_getChild(node, n), sb ); + } + } + + if (group) + { + StringBuffer_appendChar(sb, ')'); + } +} + + + +/** + * Visits the given ASTNode node. This function is really just a + * dispatcher to either SBML_formulaToString_visitFunction() or + * SBML_formulaToString_visitOther(). + */ +/* WCS addition: Change function calls to wcs modified fuctions */ +void FunctionInterfaceForJIT::FormulaFormatter_visit_wcs ( const ASTNode_t *parent, + const ASTNode_t *node, + strbuff_t *sb ) +{ + if (ASTNode_isLog10(node)) + { + StringBuffer_append(sb, "log10("); + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs(node, ASTNode_getChild(node, 1), sb); + StringBuffer_appendChar(sb, ')'); + // FormulaFormatter_visitLog10_wcs(parent, node, sb); + } + else if (ASTNode_isSqrt(node)) + { + StringBuffer_append(sb, "sqrt("); + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs(node, ASTNode_getChild(node, 1), sb); + StringBuffer_appendChar(sb, ')'); + // FormulaFormatter_visitSqrt_wcs(parent, node, sb); + } + else if (FormulaFormatter_isFunction(node)) + { + node_map::iterator nit; + LIBSBML_CPP_NAMESPACE::ASTNodeType_t node_type = ASTNode_getType(node); + nit = nodetypes.find(node_type); + if (nit != nodetypes.cend()) { + node_structure node_st = nit -> second; + ast_function_ptr node_f_ptr = std::get<0>(node_st); + if (node_f_ptr == nullptr) { + unsigned int numChildren = ASTNode_getNumChildren(node); + unsigned int n; + std::string fist_symbol = std::get<1>(node_st); + std::string middle_symbol = std::get<2>(node_st); + std::string end_symbol = std::get<3>(node_st); + StringBuffer_append(sb, fist_symbol.c_str()); + // FormulaFormatter_format(sb, node); + if (numChildren > 0) { + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs( node, LIBSBML_CPP_NAMESPACE::ASTNode_getChild(node, 0), sb ); + } + for (n = 1; n < numChildren; n++) { + StringBuffer_append(sb, middle_symbol.c_str()); + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs( node, LIBSBML_CPP_NAMESPACE::ASTNode_getChild(node, n), sb ); + } + StringBuffer_append(sb, end_symbol.c_str()); + } else { // if there is a function pointer + + } + } else { // user functions + unsigned int numChildren = ASTNode_getNumChildren(node); + unsigned int n; + StringBuffer_appendChar(sb, '('); + StringBuffer_append(sb, ASTNode_getName(node)); + + if (numChildren > 0) + { + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs( node, LIBSBML_CPP_NAMESPACE::ASTNode_getChild(node, 0), sb ); + } + + for (n = 1; n < numChildren; n++) + { + StringBuffer_appendChar(sb, ','); + StringBuffer_appendChar(sb, ' '); + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs( node, LIBSBML_CPP_NAMESPACE::ASTNode_getChild(node, n), sb ); + } + StringBuffer_appendChar(sb, ')'); + } + // FormulaFormatter_visitFunction_wcs(parent, node, sb); + } + else if (ASTNode_hasTypeAndNumChildren(node, LIBSBML_CPP_NAMESPACE::AST_MINUS, 1)) + { + StringBuffer_appendChar(sb, '-'); + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs ( node, ASTNode_getLeftChild(node), sb ); + // FormulaFormatter_visitUMinus_wcs(parent, node, sb); + } + else if (ASTNode_hasTypeAndNumChildren(node, LIBSBML_CPP_NAMESPACE::AST_PLUS, 1) || ASTNode_hasTypeAndNumChildren(node, LIBSBML_CPP_NAMESPACE::AST_TIMES, 1)) + { + FunctionInterfaceForJIT::FormulaFormatter_visit_wcs (node, ASTNode_getChild(node, 0), sb); + } + else if (ASTNode_hasTypeAndNumChildren(node, LIBSBML_CPP_NAMESPACE::AST_PLUS, 0)) + { + libsbml::StringBuffer_appendInt(sb, 0); + } + else if (ASTNode_hasTypeAndNumChildren(node, LIBSBML_CPP_NAMESPACE::AST_TIMES, 0)) + { + libsbml::StringBuffer_appendInt(sb, 1); + } + else + { + FormulaFormatter_visitOther_wcs(parent, node, sb); + } +} + +/* WCS addition: Change function call from FormulaFormatter_visit to FormulaFormatter_visit_wcs */ +std::string FunctionInterfaceForJIT::SBML_formulaToString_wcs (const ASTNode_t *tree) +{ + if (tree == nullptr) { + return nullptr; + } + + LIBSBML_CPP_NAMESPACE::StringBuffer_t *buf = LIBSBML_CPP_NAMESPACE::StringBuffer_create(1024); + + FormulaFormatter_visit_wcs(NULL, tree, buf); + std::string str(LIBSBML_CPP_NAMESPACE::StringBuffer_getBuffer(buf)); + safe_free(buf); + + return str; +} +#endif // defined(WCS_HAS_SBML) + +FunctionInterfaceForJIT::~FunctionInterfaceForJIT() +{} + +/**@}*/ +} // end of namespace wcs diff --git a/src/utils/FunctionInterfaceForJIT.hpp b/src/utils/FunctionInterfaceForJIT.hpp new file mode 100644 index 0000000..f8ca522 --- /dev/null +++ b/src/utils/FunctionInterfaceForJIT.hpp @@ -0,0 +1,74 @@ +/****************************************************************************** + * * + * 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_UTILS_FUNCTIONINTERFACEFORJIT__ +#define __WCS_UTILS_FUNCTIONINTERFACEFORJIT__ + +#if defined(WCS_HAS_CONFIG) +#include "wcs_config.hpp" +#else +#error "no config" +#endif + +#include +#include + +#if defined(WCS_HAS_SBML) +#include +#include +#endif // defined(WCS_HAS_SBML) + + +namespace wcs { +/** \addtogroup wcs_functions + * @{ */ + +#if defined(WCS_HAS_SBML) +class FunctionInterfaceForJIT { + public: + using ASTNode_t = LIBSBML_CPP_NAMESPACE::ASTNode_t; + using strbuff_t = LIBSBML_CPP_NAMESPACE::StringBuffer_t; + public: + FunctionInterfaceForJIT(); + ~FunctionInterfaceForJIT(); + + static std::string SBML_formulaToString_wcs (const ASTNode_t *tree); + + private: + static void FormulaFormatter_visit_wcs ( const ASTNode_t*parent, + const ASTNode_t *node, + strbuff_t *sb ); + + static void FormulaFormatter_visitFunction_wcs ( const ASTNode_t *parent, + const ASTNode_t *node, + strbuff_t *sb ); + + static void FormulaFormatter_visitLog10_wcs ( const ASTNode_t*parent, + const ASTNode_t *node, + strbuff_t *sb ); + + static void FormulaFormatter_visitSqrt_wcs ( const ASTNode_t *parent, + const ASTNode_t *node, + strbuff_t *sb ); + + static void FormulaFormatter_visitUMinus_wcs ( const ASTNode_t *parent, + const ASTNode_t *node, + strbuff_t *sb ); + + static void FormulaFormatter_visitOther_wcs ( const ASTNode_t *parent, + const ASTNode_t *node, + strbuff_t *sb ); + +}; +#endif // defined(WCS_HAS_SBML) + +/**@}*/ +} // end of namespace wcs +#endif // __WCS_UTILS_FUNCTIONINTERFACEFORJIT__ diff --git a/src/utils/generate_cxx_code.cpp b/src/utils/generate_cxx_code.cpp index ccb7222..e1c078c 100644 --- a/src/utils/generate_cxx_code.cpp +++ b/src/utils/generate_cxx_code.cpp @@ -31,6 +31,7 @@ #include #include #include "wcs_types.hpp" +#include "FunctionInterfaceForJIT.hpp" #include // std::thread::hardware_concurrency #if defined(WCS_HAS_SBML) @@ -170,43 +171,42 @@ generate_cxx_code::get_all_dependencies( return dependencies_no_dupl; } -static std::string add_braces_to_min( - char* math) +static double count2Concentration_converter ( + const LIBSBML_CPP_NAMESPACE::Model& model, + const LIBSBML_CPP_NAMESPACE::Species& species, + std::string& compartment, + double& comp_unit) { - std::string new_math = std::string(math); - std::string min_str ("min("); - std::string open_str ("("); - std::string close_str (")"); - std::size_t found = new_math.find(min_str); - while (found!=std::string::npos){ - new_math.insert(found, "std::"); - new_math.insert(found+9, 1 , '{'); - std::size_t close_pos = new_math.find(close_str,found+10); - std::size_t open_pos = new_math.find(open_str,found+10); - int num_open=1; - while( close_pos != std::string::npos) - { - if (open_pos == std::string::npos && num_open==1){ - new_math.insert(close_pos, 1 , '}'); - break; - } else if ( open_pos > close_pos ) { - if (num_open ==1) { - new_math.insert(close_pos, 1 , '}'); - break; - } else { - num_open = num_open - 1; - open_pos = new_math.find(open_str,close_pos+1); - close_pos = new_math.find(close_str,close_pos+1); - } - } else { - num_open = num_open + 1; - close_pos = new_math.find(close_str,open_pos+1); - open_pos = new_math.find(open_str,open_pos+1); - } - } - found = new_math.find(min_str,close_pos+1); - } - return new_math; + const LIBSBML_CPP_NAMESPACE::ListOfCompartments* compartment_list + = model.getListOfCompartments(); + const LIBSBML_CPP_NAMESPACE::ListOfUnitDefinitions* unit_definition_list + = model.getListOfUnitDefinitions(); + compartment = species.getCompartment(); + double compartment_size = compartment_list->get(species.getCompartment())->getSize(); + double avog_num = 6.02214e+23; + + + std::string compartment_unit ; + if (compartment_list->get(species.getCompartment())->isSetUnits()){ + compartment_unit = compartment_list->get(species.getCompartment())->getUnits(); + } else if (compartment_list->get(species.getCompartment())->getSpatialDimensions() == 3) { + compartment_unit = model.getVolumeUnits(); + } else if (compartment_list->get(species.getCompartment())->getSpatialDimensions() == 2) { + compartment_unit = model.getAreaUnits(); + } else if (compartment_list->get(species.getCompartment())->getSpatialDimensions() == 1) { + compartment_unit = model.getLengthUnits(); + } + const LIBSBML_CPP_NAMESPACE::UnitDefinition* unit_definition = unit_definition_list->get(compartment_unit); + const LIBSBML_CPP_NAMESPACE::ListOfUnits* unit_list + = unit_definition->getListOfUnits(); + unsigned int unitsSize = unit_list->size(); + // double comp_unit = 1.0; + comp_unit = 1.0; + for (unsigned int iu = 0u; iu < unitsSize; iu++) { + const LIBSBML_CPP_NAMESPACE::Unit* unit = unit_list->get(iu); + comp_unit = comp_unit * pow(unit->getMultiplier()*pow(10,unit->getScale()),unit->getExponent()); + } + return avog_num * compartment_size * comp_unit; } @@ -228,7 +228,7 @@ generate_cxx_code::return_denominators( } else { if ( std::find(denominators.cbegin(), denominators.cend(), SBML_formulaToString(denominator)) == denominators.cend()) { - std::string new_denominator = add_braces_to_min(SBML_formulaToString(denominator)); + std::string new_denominator = FunctionInterfaceForJIT::SBML_formulaToString_wcs(denominator); denominators.push_back(new_denominator); } } @@ -406,8 +406,15 @@ void generate_cxx_code::find_used_params( // Find used parameters in the rates for (unsigned int ic = 0u; ic < num_reactions; ic++) { const LIBSBML_CPP_NAMESPACE::Reaction& reaction = *(reaction_list->get(ic)); + // const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list + // = reaction.getKineticLaw()->getListOfLocalParameters(); + // unsigned int num_localparameters = local_parameter_list->size(); + // //genfile << "reaction" << reaction.getIdAttribute() <<": "; using reaction_parameters_t = std::unordered_set; reaction_parameters_t lpset; + // for (unsigned int pi = 0u; pi < num_localparameters; pi++) { + // const LIBSBML_CPP_NAMESPACE::LocalParameter* localparameter = local_parameter_list->get(pi); + // lpset.insert(localparameter->getIdAttribute()); // Check SBML level for local parameters if (model.getLevel() > 2) { const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list @@ -1123,6 +1130,13 @@ void generate_cxx_code::print_reaction_rates( genfile << "#include \"" + header + '"' + "\n\n"; for (unsigned int ic = rid_start; ic < rid_end; ic++) { const LIBSBML_CPP_NAMESPACE::Reaction& reaction = *(reaction_list->get(ic)); + // const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list + // = reaction.getKineticLaw()->getListOfLocalParameters(); + // unsigned int num_localparameters = local_parameter_list->size(); + // a map to keep the compartment and their units + using compartment_t = std::unordered_map ; + typename compartment_t::const_iterator rci; + compartment_t reaction_comp; unsigned int num_localparameters = 0u; if (model.getLevel() > 2) { const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list @@ -1134,6 +1148,7 @@ void generate_cxx_code::print_reaction_rates( num_localparameters = local_parameter_list->size(); } + genfile << "extern \"C\" " << Real << " wcs__rate_" << reaction.getIdAttribute() << "(const std::vector<" << Real <<">& __input) {\n"; @@ -1141,13 +1156,18 @@ void generate_cxx_code::print_reaction_rates( //print reaction's local parameters using reaction_local_parameters_t = std::unordered_set; reaction_local_parameters_t localpset; + // for (unsigned int pi = 0u; pi < num_localparameters; pi++) { + // const LIBSBML_CPP_NAMESPACE::LocalParameter* localparameter = local_parameter_list->get(pi); + // genfile << " constexpr " << Real << " " << localparameter->getIdAttribute() << " = " + // << localparameter->getValue() << ";\n"; + // localpset.insert(localparameter->getIdAttribute()); if (model.getLevel() > 2) { const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list = reaction.getKineticLaw()->getListOfLocalParameters(); for (unsigned int pi = 0u; pi < num_localparameters; pi++) { const LIBSBML_CPP_NAMESPACE::LocalParameter* localparameter = local_parameter_list->get(pi); - //genfile << " constexpr " << Real << " " << localparameter->getIdAttribute() << " = " - // << localparameter->getValue() << ";\n"; + genfile << " constexpr " << Real << " " << localparameter->getIdAttribute() << " = " + << localparameter->getValue() << ";\n"; localpset.insert(localparameter->getIdAttribute()); } } else { @@ -1155,8 +1175,8 @@ void generate_cxx_code::print_reaction_rates( = reaction.getKineticLaw()->getListOfParameters(); for (unsigned int pi = 0u; pi < num_localparameters; pi++) { const LIBSBML_CPP_NAMESPACE::Parameter* localparameter = local_parameter_list->get(pi); - //genfile << " constexpr " << Real << " " << localparameter->getIdAttribute() << " = " - // << localparameter->getValue() << ";\n"; + genfile << " constexpr " << Real << " " << localparameter->getIdAttribute() << " = " + << localparameter->getValue() << ";\n"; localpset.insert(localparameter->getIdAttribute()); } } @@ -1192,12 +1212,15 @@ void generate_cxx_code::print_reaction_rates( for (const auto& x: input_map) { var_names_ord[x.second] = x.first; } - //print first parameters in formula taking input std::vector::const_iterator it; std::vector par_names, par_names_nf; all_var_names = var_names; int par_index = 0; + const ListOfSpecies* species_list = model.getListOfSpecies(); + double avog_num = 6.02214e+23; + bool concentration_used = false; + for (it = var_names_ord.cbegin(); it < var_names_ord.cend(); it++){ rrdit = rate_rules_dep_map.find(*it); evassigit = ev_assign.find(*it); @@ -1215,10 +1238,38 @@ void generate_cxx_code::print_reaction_rates( all_var_names_it = all_var_names.find(*itf); if (all_var_names_it == all_var_names.cend()) { genfile << " " << Real << " " << *itf << " = __input[" << par_index++ << "];\n"; + if (species_list->get(*itf) != NULL) { + if (species_list->get(*itf)->isSetInitialConcentration()) { + double comp_unit; + std::string compart_name; + double convert_factor = count2Concentration_converter(model, *species_list->get(*itf), compart_name, comp_unit); + // Add values for the compartment in map + rci = reaction_comp.find(compart_name) ; + if (rci == reaction_comp.end()) { + reaction_comp.insert(std::make_pair(compart_name, comp_unit)); + } + genfile << " " << *itf << " = " << *itf << " / " << convert_factor << ";\n"; + concentration_used = true; + } + } var_names.erase(*itf); } else { if (var_names_it != var_names.cend()) { genfile << " " << Real << " " << *itf << " = __input[" << par_index++ << "];\n"; + if (species_list->get(*itf) != NULL) { + if (species_list->get(*itf)->isSetInitialConcentration()) { + double comp_unit; + std::string compart_name; + double convert_factor = count2Concentration_converter(model, *species_list->get(*itf), compart_name, comp_unit); + // Add values for the compartment in map + rci = reaction_comp.find(compart_name) ; + if (rci == reaction_comp.end()) { + reaction_comp.insert(std::make_pair(compart_name, comp_unit)); + } + genfile << " " << *itf << " = " << *itf << " / " << convert_factor << ";\n"; + concentration_used = true; + } + } var_names.erase(*itf); } } @@ -1233,7 +1284,23 @@ void generate_cxx_code::print_reaction_rates( genfile << " wcs_global_var." << *it << " = wcs__rate_" << *it << "(" << function_input << ");\n"; } else { genfile << " " << Real << " " << *it << " = __input[" << par_index++ << "];\n"; - sum_stoich = sum_stoich + reaction.getReactant(*it)->getStoichiometry(); + if (species_list->get(*it) != NULL) { + if (species_list->get(*it)->isSetInitialConcentration()) { + double comp_unit; + std::string compart_name; + double convert_factor = count2Concentration_converter(model, *species_list->get(*it), compart_name, comp_unit); + // Add values for the compartment in map + rci = reaction_comp.find(compart_name) ; + if (rci == reaction_comp.end()) { + reaction_comp.insert(std::make_pair(compart_name, comp_unit)); + } + genfile << " " << *it << " = " << *it << " / " << convert_factor << ";\n"; + concentration_used = true; + } + } + if (reaction.getReactant(*it) != NULL) { + sum_stoich = sum_stoich + reaction.getReactant(*it)->getStoichiometry(); + } par_names.push_back(*it); } var_names.erase(*it); @@ -1276,172 +1343,6 @@ void generate_cxx_code::print_reaction_rates( par_names_nf.push_back(x); } } - const LIBSBML_CPP_NAMESPACE::ListOfUnitDefinitions* unit_definition_list - = model.getListOfUnitDefinitions(); - if (model.getLevel() > 2) { - int has_only_substance = 0; - double compartment_size = 0; - const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list - = reaction.getKineticLaw()->getListOfLocalParameters(); - const LIBSBML_CPP_NAMESPACE::ListOfCompartments* compartment_list - = model.getListOfCompartments(); - for (unsigned int pi = 0u; pi < num_localparameters; pi++) { - const LIBSBML_CPP_NAMESPACE::LocalParameter* localparameter = local_parameter_list->get(pi); - - // Check hasOnlySubstance - const ListOfSpecies* species_list = model.getListOfSpecies(); - const unsigned int num_species = species_list->size(); - for (unsigned int ic = 0u; ic < num_species; ic++) { - if (species_list->get(ic)->getHasOnlySubstanceUnits()) { - has_only_substance = 1; - } else { - compartment_size = compartment_list->get(species_list->get(ic)->getCompartment())->getSize(); - } - } - - if (localparameter->isSetUnits()){ - std::string local_par_unit = localparameter->getUnits(); - const LIBSBML_CPP_NAMESPACE::UnitDefinition* unit_definition = unit_definition_list->get(local_par_unit); - const LIBSBML_CPP_NAMESPACE::ListOfUnits* unit_list - = unit_definition->getListOfUnits(); - unsigned int unitsSize = unit_list->size(); - double comp_unit = 1.0; - for (unsigned int iu = 0u; iu < unitsSize; iu++) { - const LIBSBML_CPP_NAMESPACE::Unit* unit = unit_list->get(iu); - comp_unit = comp_unit * pow(unit->getMultiplier()*pow(10,unit->getScale()),unit->getExponent()); - } - genfile << " constexpr " << Real << " " << localparameter->getIdAttribute() << " = " - << localparameter->getValue() * comp_unit << ";\n"; - } else { - //Volume units/(substance units *timeunits) - if (model.isSetVolumeUnits() && model.isSetSubstanceUnits() && model.isSetTimeUnits()){ - std::string volume_units = model.getVolumeUnits(); - std::string substance_units = model.getSubstanceUnits(); - std::string time_units = model.getTimeUnits(); - double comp_unit = 1.0; - const LIBSBML_CPP_NAMESPACE::UnitDefinition* unit_definition; - const LIBSBML_CPP_NAMESPACE::ListOfUnits* unit_list; - unsigned int unitsSize = 0u; - - if (unit_definition_list->get(volume_units) != NULL) { - unit_definition = unit_definition_list->get(volume_units); - unit_list = unit_definition->getListOfUnits(); - unitsSize = unit_list->size(); - for (unsigned int iu = 0u; iu < unitsSize; iu++) { - const LIBSBML_CPP_NAMESPACE::Unit* unit = unit_list->get(iu); - comp_unit = comp_unit * pow(unit->getMultiplier()*pow(10,unit->getScale()),(sum_stoich-1) * unit->getExponent()); - } - } - - if (unit_definition_list->get(substance_units) != NULL) { - unit_definition = unit_definition_list->get(substance_units); - unit_list = unit_definition->getListOfUnits(); - unitsSize = unit_list->size(); - for (unsigned int iu = 0u; iu < unitsSize; iu++) { - const LIBSBML_CPP_NAMESPACE::Unit* unit = unit_list->get(iu); - if (unit->getScale() == -3) { - comp_unit = comp_unit * pow(unit->getMultiplier(),-((sum_stoich-1) * unit->getExponent())); - } else { - comp_unit = comp_unit * pow(unit->getMultiplier()*pow(10,unit->getScale()),-((sum_stoich-1) * unit->getExponent())); - } - } - } - - if (unit_definition_list->get(time_units) != NULL) { - unit_definition = unit_definition_list->get(time_units); - unit_list = unit_definition->getListOfUnits(); - unitsSize = unit_list->size(); - for (unsigned int iu = 0u; iu < unitsSize; iu++) { - const LIBSBML_CPP_NAMESPACE::Unit* unit = unit_list->get(iu); - comp_unit = comp_unit * pow(unit->getMultiplier()*pow(10,unit->getScale()),-((sum_stoich-1) * unit->getExponent())); - } - } - if (has_only_substance ==1) { - // genfile << " constexpr " << Real << " " << localparameter->getIdAttribute() << " = " - // << localparameter->getValue() * comp_unit << ";\n"; - genfile << " constexpr " << Real << " " << localparameter->getIdAttribute() << " = " - << localparameter->getValue() << ";\n"; - } else { - // genfile << " constexpr " << Real << " " << localparameter->getIdAttribute() << " = " - // << localparameter->getValue() * comp_unit / compartment_size << ";\n"; - // genfile << " constexpr " << Real << " " << localparameter->getIdAttribute() << " = " - // << localparameter->getValue() * comp_unit << ";\n"; - genfile << " constexpr " << Real << " " << localparameter->getIdAttribute() << " = " - << localparameter->getValue() << ";\n"; - } - } - } - } - } else { - const LIBSBML_CPP_NAMESPACE::ListOfParameters* local_parameter_list - = reaction.getKineticLaw()->getListOfParameters(); - for (unsigned int pi = 0u; pi < num_localparameters; pi++) { - const LIBSBML_CPP_NAMESPACE::Parameter* localparameter = local_parameter_list->get(pi); - - if (localparameter->isSetUnits()){ - std::string local_par_unit = localparameter->getUnits(); - const LIBSBML_CPP_NAMESPACE::UnitDefinition* unit_definition = unit_definition_list->get(local_par_unit); - const LIBSBML_CPP_NAMESPACE::ListOfUnits* unit_list - = unit_definition->getListOfUnits(); - unsigned int unitsSize = unit_list->size(); - double comp_unit = 1.0; - for (unsigned int iu = 0u; iu < unitsSize; iu++) { - const LIBSBML_CPP_NAMESPACE::Unit* unit = unit_list->get(iu); - comp_unit = comp_unit * pow(unit->getMultiplier()*pow(10,unit->getScale()),unit->getExponent()); - } - genfile << " constexpr " << Real << " " << localparameter->getIdAttribute() << " = " - << localparameter->getValue() * comp_unit << ";\n"; - } else { - //Volume units/(substance units *timeunits) - if (model.isSetVolumeUnits() && model.isSetSubstanceUnits() && model.isSetTimeUnits()){ - std::string volume_units = model.getVolumeUnits(); - std::string substance_units = model.getSubstanceUnits(); - std::string time_units = model.getTimeUnits(); - double comp_unit = 1.0; - const LIBSBML_CPP_NAMESPACE::UnitDefinition* unit_definition; - const LIBSBML_CPP_NAMESPACE::ListOfUnits* unit_list; - unsigned int unitsSize = 0u; - - if (unit_definition_list->get(volume_units) != NULL) { - unit_definition = unit_definition_list->get(volume_units); - unit_list = unit_definition->getListOfUnits(); - unitsSize = unit_list->size(); - for (unsigned int iu = 0u; iu < unitsSize; iu++) { - const LIBSBML_CPP_NAMESPACE::Unit* unit = unit_list->get(iu); - comp_unit = comp_unit * pow(unit->getMultiplier()*pow(10,unit->getScale()),(sum_stoich-1) * unit->getExponent()); - } - } - - if (unit_definition_list->get(substance_units) != NULL) { - unit_definition = unit_definition_list->get(substance_units); - unit_list = unit_definition->getListOfUnits(); - unitsSize = unit_list->size(); - for (unsigned int iu = 0u; iu < unitsSize; iu++) { - const LIBSBML_CPP_NAMESPACE::Unit* unit = unit_list->get(iu); - if (unit->getScale() == -3) { - comp_unit = comp_unit * pow(unit->getMultiplier(),-((sum_stoich-1) * unit->getExponent())); - } else { - comp_unit = comp_unit * pow(unit->getMultiplier()*pow(10,unit->getScale()),-((sum_stoich-1) * unit->getExponent())); - } - } - } - - if (unit_definition_list->get(time_units) != NULL) { - unit_definition = unit_definition_list->get(time_units); - unit_list = unit_definition->getListOfUnits(); - unitsSize = unit_list->size(); - for (unsigned int iu = 0u; iu < unitsSize; iu++) { - const LIBSBML_CPP_NAMESPACE::Unit* unit = unit_list->get(iu); - comp_unit = comp_unit * pow(unit->getMultiplier()*pow(10,unit->getScale()),-((sum_stoich-1) * unit->getExponent())); - } - } - - genfile << " constexpr " << Real << " " << localparameter->getIdAttribute() << " = " - << localparameter->getValue() * comp_unit << ";\n"; - } - } - } - } // put input parameters into maps dep_params_f.insert(std::make_pair(reaction.getIdAttribute(), par_names)); @@ -1474,7 +1375,7 @@ void generate_cxx_code::print_reaction_rates( } math = *reaction.getKineticLaw()->getMath(); update_scope_ast_node(math, wcs_all_var, wcs_all_const, localpset); - std::string new_math = add_braces_to_min(SBML_formulaToString(&math)); + std::string new_math = FunctionInterfaceForJIT::SBML_formulaToString_wcs(&math); genfile << " " << Real << " " << reaction.getIdAttribute() << " = " << new_math //SBML_formulaToString( &math) << ";\n"; @@ -1532,7 +1433,16 @@ void generate_cxx_code::print_reaction_rates( genfile << " }\n"; //genfile << " printf(\" Result of %s : %f \", \"" << reaction.getIdAttribute() << "\"," // << reaction.getIdAttribute() << ");\n"; - genfile << " return " << reaction.getIdAttribute() << ";\n"; + if (!concentration_used) { + genfile << " return " << reaction.getIdAttribute() << ";\n"; + } else { + // Calculate the compartment unit for all the compartments + double comp_unit = 1.0; + for ( auto rci = reaction_comp.begin(); rci != reaction_comp.end(); ++rci ){ + comp_unit = comp_unit * rci->second; + } + genfile << " return " << reaction.getIdAttribute() << " * " << avog_num << " * " << comp_unit << ";\n"; + } genfile << "}\n\n"; } } From 8e8d0bac8ddcfdc036aa788aae30e30ca5029fab Mon Sep 17 00:00:00 2001 From: Konstantia Georgouli Date: Tue, 24 Aug 2021 10:39:20 -0700 Subject: [PATCH 03/14] Support hybrid simulation mode with ROSS --- CMakeLists.txt | 19 + src/CMakeLists.txt | 2 + src/hybrid.cpp | 444 ++++++++++++++++++ src/hybrid.hpp | 53 +++ src/ross_files/CMakeLists.txt | 36 ++ src/ross_files/Checkpoint.hh | 46 ++ src/ross_files/Funnel.hh | 560 +++++++++++++++++++++++ src/ross_files/actor.cc | 13 + src/ross_files/actor.hh | 138 ++++++ src/ross_files/generateInfrastructure.py | 188 ++++++++ src/ross_files/simtime.hh | 10 + 11 files changed, 1509 insertions(+) create mode 100644 src/hybrid.cpp create mode 100644 src/hybrid.hpp create mode 100644 src/ross_files/CMakeLists.txt create mode 100644 src/ross_files/Checkpoint.hh create mode 100644 src/ross_files/Funnel.hh create mode 100644 src/ross_files/actor.cc create mode 100644 src/ross_files/actor.hh create mode 100644 src/ross_files/generateInfrastructure.py create mode 100644 src/ross_files/simtime.hh diff --git a/CMakeLists.txt b/CMakeLists.txt index 8e3c745..899f545 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -620,6 +620,19 @@ if (WCS_HAS_ROSS) list(APPEND WCS_EXEC_TARGETS des-bin) endif (WCS_HAS_ROSS) +# add executable hybrid +if (WCS_HAS_ROSS) + add_executable( hybrid-bin src/hybrid.cpp ) + target_include_directories(hybrid-bin SYSTEM PUBLIC ${ROSS_INCLUDE_DIRS}) + target_include_directories(hybrid-bin PUBLIC + $ + $ + $) + target_link_libraries(hybrid-bin PRIVATE ROSS::ROSS MPI::MPI_CXX wcs ${LIB_FILESYSTEM}) + set_target_properties(hybrid-bin PROPERTIES OUTPUT_NAME hybrid) + list(APPEND WCS_EXEC_TARGETS hybrid-bin) +endif (WCS_HAS_ROSS) + # Build tests add_executable( t_state_rngen-bin src/utils/unit_tests/t_state_rngen.cpp ) target_include_directories(t_state_rngen-bin PUBLIC @@ -713,6 +726,7 @@ set(INCLUDE_INSTALL_DIRS "${CMAKE_SOURCE_DIR}/src/proto" "${CMAKE_SOURCE_DIR}/src/reaction_network" "${CMAKE_SOURCE_DIR}/src/sim_methods" + "${CMAKE_SOURCE_DIR}/src/ross_files" "${CMAKE_SOURCE_DIR}/src/utils") set(LIB_INSTALL_DIR "${CMAKE_BINARY_DIR}") set(EXTRA_CMAKE_MODULE_DIR "${CMAKE_SOURCE_DIR}/cmake/modules") @@ -787,11 +801,16 @@ install( DIRECTORY "${PROJECT_SOURCE_DIR}/src/utils" DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}" FILES_MATCHING PATTERN "*.h" PATTERN "*.hpp") +install( + DIRECTORY "${PROJECT_SOURCE_DIR}/src/ross_files" + DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}" + FILES_MATCHING PATTERN "*.h" PATTERN "*.hpp" PATTERN "*.hh") install( FILES "${PROJECT_BINARY_DIR}/wcs_config.hpp" "${PROJECT_SOURCE_DIR}/src/wcs_types.hpp" "${PROJECT_SOURCE_DIR}/src/bgl.hpp" "${PROJECT_SOURCE_DIR}/src/des.hpp" + "${PROJECT_SOURCE_DIR}/src/hybrid.hpp" "${PROJECT_SOURCE_DIR}/src/wcs-ross-bf.hpp" DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}") diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index c55fc58..ab109fc 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -3,6 +3,7 @@ set_full_path(THIS_DIR_HEADERS bgl.hpp wcs_types.hpp des.hpp + hybrid.hpp wcs-ross-bf.hpp ) @@ -16,6 +17,7 @@ add_subdirectory(proto) add_subdirectory(reaction_network) add_subdirectory(sim_methods) add_subdirectory(utils) +add_subdirectory(ross_files) # Propagate the files up the tree set(WCS_HEADERS "${WCS_HEADERS}" "${THIS_DIR_HEADERS}" PARENT_SCOPE) diff --git a/src/hybrid.cpp b/src/hybrid.cpp new file mode 100644 index 0000000..f67582c --- /dev/null +++ b/src/hybrid.cpp @@ -0,0 +1,444 @@ +/****************************************************************************** + * * + * 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 * + * * + ******************************************************************************/ + +#if defined(WCS_HAS_CONFIG) +#include "wcs_config.hpp" +#else +#error "no config" +#endif + +#if defined(WCS_HAS_ROSS) + +#include +#include // memset +#include +#include +#include +#include "utils/write_graphviz.hpp" +#include "utils/timer.hpp" +#include "utils/to_string.hpp" +#include "reaction_network/network.hpp" +#include "hybrid.hpp" +#include "ross_files/simtime.hh" +#include "ross_files/Funnel.hh" + +#define OPTIONS "hi:s:t:" +static const struct option longopts[] = { + {"help", no_argument, 0, 'h'}, + {"interval", required_argument, 0, 'i'}, + {"seedbase", required_argument, 0, 's'}, + {"trust_region", required_argument, 0, 't'}, + { 0, 0, 0, 0 }, +}; + +void do_nothing(void *,void *,void *) {} + +tw_peid ctr_map(tw_lpid gid) { + return (tw_peid) gid / g_tw_nlp; +} + +template +class RegisterBufferAndOffset { + public: + static std::vector buffer; + static int offset; + + static void init_f(void *state, tw_lp* me) + { + LP::init_lp(state, me); + buffer[me->gid-offset] = reinterpret_cast(state); + } + + static void init(const int bufferSize, const int _offset) + { + buffer.resize(bufferSize); + offset = _offset; + } + static TTT** getBuffer() { + return &buffer[0]; + } +}; + +tw_lptype Funnellps[] = { + { (init_f) RegisterBufferAndOffset::init_f, + (pre_run_f) do_nothing, + (event_f) Funnel::forward_event, + (revent_f) Funnel::backward_event, + (commit_f) Funnel::commit_event, + (final_f) Funnel::final_lp, + (map_f) ctr_map, + //NULL, + sizeof(Funnel) + }, + {0} +}; +template<> +std::vector RegisterBufferAndOffset::buffer = std::vector(); +template<> +int RegisterBufferAndOffset::offset = 0; + +tw_lptype NRMlps[] = { + { (init_f) RegisterBufferAndOffset::init_f, + (pre_run_f) do_nothing, + (event_f) NextReactionMethodCrucible::forward_event, + (revent_f) NextReactionMethodCrucible::backward_event, + (commit_f) NextReactionMethodCrucible::commit_event, + (final_f) NextReactionMethodCrucible::final_lp, + (map_f) ctr_map, + //NULL, + sizeof(NextReactionMethodCrucible) + }, + {0} +}; + +template<> +std::vector RegisterBufferAndOffset::buffer = std::vector(); +template<> +int RegisterBufferAndOffset::offset = 0; + + +tw_lptype Heartbeatlps[] = { + { (init_f) RegisterBufferAndOffset::init_f, + (pre_run_f) do_nothing, + (event_f) Heartbeat::forward_event, + (revent_f) Heartbeat::backward_event, + (commit_f) Heartbeat::commit_event, + (final_f) Heartbeat::final_lp, + (map_f) ctr_map, + //NULL, + sizeof(Heartbeat) + }, + {0} +}; + + +template<> +std::vector RegisterBufferAndOffset::buffer = std::vector(); +template<> +int RegisterBufferAndOffset::offset = 0; + +// void post_lpinit_setup(tw_pe* pe); + +tw_petype MyPEType = { + NULL, + post_lpinit_setup, + NULL, + NULL +}; + + + +//FIXME, use graph ids +enum ReactionEnum { + A_TO_X, + Y_2X_TO_3X, + B_X_TO_Y_D, + X_TO_E, + numReactions +}; + +enum SpeciesEnum { + A_idx, + B_idx, + D_idx, + E_idx, + X_idx, + Y_idx, + numSpecies +}; + + + +// FIXME, comment out, every stoic reference has to check the graph +std::vector> stoic(numReactions); + + +//FIXME, all of these functions are linked in with dlsym trick + +int factor = 200; +Real k_A_TO_X = 1; +Real k_Y_2X_TO_3X = 2.7574E-06/factor/factor; +Real k_B_X_TO_Y_D = 0.001660538/factor; +Real k_X_TO_E = 1; + +Real function_A_TO_X(const std::vector& species) { + // for (auto specie : species){ + // std::cout << specie << std::endl; + // } + return k_A_TO_X*species[0]; + +} + +Real function_Y_2X_TO_3X(const std::vector& species) { + return k_Y_2X_TO_3X*species[0]*species[0]*species[1]; +} + +Real function_B_X_TO_Y_D(const std::vector& species) { + return k_B_X_TO_Y_D*species[0]*species[1]; +} + +Real function_X_TO_E(const std::vector& species) { + return k_X_TO_E*species[0]; +} + + +void print_usage(const std::string exec, int code) +{ + std::cerr << + "Usage: " << exec << " .xml\n" + " Run the stochastic simulation algorithm (SSA) on a reaction\n" + " network graph (.xml).\n" + "\n" + " OPTIONS:\n" + " -h, --help\n" + " Display this usage information\n" + "\n" + " -i, --interval\n" + " Specify the interval\n" + "\n" + " -s, --seedbase\n" + " Specify the seed for random number generator\n" + "\n" + " -t, --trust_region\n" + " Specify the trust region\n" + "\n"; + exit(code); +} + + +int main(int argc, char **argv) +{ + // get rid of error if compiled w/ MEMORY queues + //g_tw_memory_nqueues=1; + //g_tw_memory_nqueues = 16; // give at least 16 memory queue event + + //printf("Calling tw_opt_add\n"); + //tw_opt_add(app_opt); + + + int c; + int rc = EXIT_SUCCESS; + std::string outfile; + double interval = 0.2; + int seedBase = 137; + int trust_region = 20000; + while ((c = getopt_long(argc, argv, OPTIONS, longopts, NULL)) != -1) { + switch (c) { + case 'h': /* --help */ + print_usage(argv[0], 0); + break; + case 'i': /* --interval */ + interval = atof(optarg); + break; + case 's': /* --seedbase */ + seedBase = static_cast(atoi(optarg)); + break; + case 't': /* --trust_region */ + trust_region = static_cast(atoi(optarg)); + break; + default: + print_usage(argv[0], 1); + break; + } + } + if (optind != (argc - 1)) { + print_usage (argv[0], 1); + } + std::string fn(argv[optind]); + std::shared_ptr rnet_ptr = std::make_shared(); + wcs::Network& rnet = *rnet_ptr; + rnet.load(fn); + rnet.init(); + const wcs::Network::graph_t& g = rnet.graph(); + + + //FIXME, needs to come from graph => DONE + int numReactions = rnet.reaction_list().size(); + + //printf("Calling tw_init\n"); + tw_init(&argc, &argv); + //printf("tw_init returned.\n"); + // g_tw_nlp = numReactions+numReactions+numReactions; + g_tw_nlp = numReactions+numReactions+1; + g_tw_events_per_pe = g_tw_nlp * 100; + + //FIXME, yeom2 . I need to fill out these arrays with pointers to the actual objects. + RegisterBufferAndOffset::init(numReactions,0); + RegisterBufferAndOffset::init(numReactions,numReactions); + RegisterBufferAndOffset::init(1,2*numReactions); + + tw_define_lps(g_tw_nlp, messageSize()); + + for(int ii = 0; ii < numReactions ; ii++) { + tw_lp_settype(ii, &Funnellps[0]); + } + for(int ii = numReactions; ii < numReactions+numReactions; ii++) { + tw_lp_settype(ii, &NRMlps[0]); + } + tw_lp_settype(2*numReactions, &Heartbeatlps[0]); + + tw_pe_settype(&MyPEType); + // std::cout << "sampling end: " << g_st_sampling_end < 0) { + g_tw_ts_end = g_st_sampling_end; // default 100000 + } else { + g_tw_ts_end = 1; // default 100000 + } + + // g_tw_events_per_pe =10; //default 800 + if( g_tw_mynode == 0 ) + { + std::cout << "=========================================" << std::endl; + std::cout << "WCS ROSS Configuration.............." << std::endl; + // std::cout << " run_id:\t" + std::string(run_id) << std::endl; + std::cout << " nlp_per_pe:\t" << g_tw_nlp << std::endl; + std::cout << " g_tw_ts_end:\t" << g_tw_ts_end << std::endl; + std::cout << " g_tw_events_per_pe:\t" << g_tw_events_per_pe << std::endl; + + std::cout << " gvt-interval:\t" << g_tw_gvt_interval << std::endl;; + std::cout << " extramem:\t" << g_tw_events_per_pe_extra << std::endl; + std::cout << " ......................................" << std::endl; + std::cout << " Num nodes:\t" << tw_nnodes() << std::endl; + // std::cout << " Message size:\t" << sizeof(WCS_Message) << std::endl; + std::cout << "=========================================" + << std::endl << std::endl; + } + + + tw_run(); + + tw_end(); + +// return 0; +return EXIT_SUCCESS; +} + +void post_lpinit_setup(tw_pe* pe) { + + Funnel** funnels = RegisterBufferAndOffset::getBuffer(); + NextReactionMethodCrucible** crucibles = RegisterBufferAndOffset::getBuffer(); + Heartbeat** heartbeats = RegisterBufferAndOffset::getBuffer(); + + + //FIXME, output interval comes from options? + double interval = 0.2; + heartbeats[0]->setInterval(interval); + + //FIXME, seed comes from options + int seedBase = 137; + for (int ireaction=0; ireactionsetTag(ireaction); + crucibles[ireaction]->setSeed(seedBase+ireaction); + + funnels[ireaction]->setTag(ireaction); + funnels[ireaction]->setSeed(seedBase+numReactions+ireaction); + } + + + //FIXME, inititial values are looked up in graph + //setup initial conditions + SpeciesValue initial_value[numSpecies]; + initial_value[A_idx] = 301*factor; + initial_value[B_idx] = 1806*factor; + initial_value[D_idx] = 0*factor; + initial_value[E_idx] = 0*factor; + initial_value[X_idx] = 1806*factor; + initial_value[Y_idx] = 1806*factor; + + //FIXME, numerical parameter, comes from options + int trust_region = 20000; + + + //FIXME, run addSpecies for each dep species, add the rate function. + + //stoic[A_TO_X][A_idx] = -1; + stoic[A_TO_X][X_idx] = 1; + funnels[A_TO_X]->_rateFunction = &function_A_TO_X; + funnels[A_TO_X]->addSpecies(A_idx, initial_value[A_idx], trust_region); + //funnels[A_TO_X]->addSpecies(X_idx, initial_value[X_idx]); + + stoic[Y_2X_TO_3X][X_idx] = -2+3; + stoic[Y_2X_TO_3X][Y_idx] = -1; + funnels[Y_2X_TO_3X]->_rateFunction = function_Y_2X_TO_3X; + funnels[Y_2X_TO_3X]->addSpecies(X_idx, initial_value[X_idx], trust_region); + funnels[Y_2X_TO_3X]->addSpecies(Y_idx, initial_value[Y_idx], trust_region); + + //stoic[B_X_TO_Y_D][B_idx] = -1; + stoic[B_X_TO_Y_D][X_idx] = -1; + stoic[B_X_TO_Y_D][Y_idx] = 1; + //stoic[B_X_TO_Y_D][D_idx] = 1; + funnels[B_X_TO_Y_D]->_rateFunction = function_B_X_TO_Y_D; + funnels[B_X_TO_Y_D]->addSpecies(B_idx, initial_value[B_idx], trust_region); + //funnels[B_X_TO_Y_D]->addSpecies(D_idx, initial_value[D_idx], trust_region); + funnels[B_X_TO_Y_D]->addSpecies(X_idx, initial_value[X_idx], trust_region); + //funnels[B_X_TO_Y_D]->addSpecies(Y_idx, initial_value[Y_idx], trust_region); + + stoic[X_TO_E][X_idx] = -1; + //stoic[X_TO_E][E_idx] = 1; + funnels[X_TO_E]->_rateFunction = function_X_TO_E; + //funnels[X_TO_E]->addSpecies(E_idx, initial_value[E_idx], trust_region); + funnels[X_TO_E]->addSpecies(X_idx, initial_value[X_idx], trust_region); + + //Connect up the funnels to the crucibles they own. + for( int ii=0; ii speciesForThisRxn; + //replace stoic instance, possibly do better than a n^2 loop because we have a graph + for (auto iter : stoic[ireaction]) { + speciesForThisRxn.push_back(iter.first); + } + funnels[ireaction]->setupRecvFromReaction(ireaction); + for (int jreaction=0; jreactiondependsOnSpecies(speciesTag)) { + CONNECT(*(crucibles[ireaction]),NextReactionMethodCrucible::fire, + *(funnels[jreaction]),Funnel::firedReaction); + if (ireaction != jreaction) { + funnels[ireaction]->setupSendToReaction(jreaction, funnels[jreaction]->getID()); + funnels[jreaction]->setupRecvFromReaction(ireaction); + } + break; + } + } + } + } + + //connect the printing to the heartbeat for output + CONNECT((*heartbeats[0]),Heartbeat::activate,(*heartbeats[0]),Heartbeat::activated); + for (int ireaction=0; ireactionbeginReactions(tnow); + + } + NoopMsg msg; + heartbeats[0]->activate(tnow, msg); +} + + + + +#endif // defined(WCS_HAS_ROSS) diff --git a/src/hybrid.hpp b/src/hybrid.hpp new file mode 100644 index 0000000..ec12ac3 --- /dev/null +++ b/src/hybrid.hpp @@ -0,0 +1,53 @@ +/****************************************************************************** + * * + * 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_HYBRID_HPP +#define WCS_HYBRID_HPP + +#if defined(WCS_HAS_CONFIG) +#include "wcs_config.hpp" +#else +#error "no config" +#endif + +#if defined(WCS_HAS_ROSS) + +#include +#include +#include +#include "reaction_network/network.hpp" +#include "sim_methods/ssa_nrm.hpp" +#include "params/ssa_params.hpp" + + + +void do_nothing(void *,void *,void *); +tw_peid ctr_map(tw_lpid gid); +void post_lpinit_setup(tw_pe* pe); + + + + +// //----------------- +// // ROSS options +// //----------------- +// static unsigned int nlp_per_pe = 1u; +// static char run_id[1024] = "wcs-ross"; + +// const tw_optdef app_opt[] = +// { +// TWOPT_GROUP("WCS ROSS"), +// TWOPT_UINT("nlp", nlp_per_pe, "Number of LPs per processor"), +// TWOPT_CHAR("run", run_id, "User supplied run name"), +// TWOPT_END() +// }; + +#endif // defined(WCS_HAS_ROSS) +#endif // WCS_HYBRID_HPP diff --git a/src/ross_files/CMakeLists.txt b/src/ross_files/CMakeLists.txt new file mode 100644 index 0000000..e8c70ed --- /dev/null +++ b/src/ross_files/CMakeLists.txt @@ -0,0 +1,36 @@ +# blt_add_executable( +# NAME ross_files +# SOURCES actor.cc +# DEPENDS_ON ROSS mpi +# ) +# find_package(Python COMPONENTS Interpreter) + +ADD_CUSTOM_COMMAND(OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/actor.gen.cc + COMMAND ${Python_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/generateInfrastructure.py ${actor_srcs} > ${CMAKE_CURRENT_BINARY_DIR}/actor.gen.cc + MAIN_DEPENDENCY generateInfrastructure.py + DEPENDS generateInfrastructure.py ${actor_srcs} + VERBATIM) +SET_PROPERTY(SOURCE ${CMAKE_CURRENT_BINARY_DIR}/actor.gen.cc + PROPERTY INCLUDE_DIRECTORIES ${CMAKE_CURRENT_SOURCE_DIR}) + + + +# Add the header and source files for this directory +# set_full_path(THIS_DIR_HEADERS +# actor.hh +# Checkpoint.hh +# Funnel.hh +# simtime.hh +# ) + +# set_full_path(THIS_DIR_SOURCES +# actor.cc +# ) + +# # Add the subdirectories + +# # Propagate the files up the tree +# set(WCS_HEADERS "${WCS_HEADERS}" "${THIS_DIR_HEADERS}" PARENT_SCOPE) +# set(WCS_SOURCES "${WCS_SOURCES}" "${THIS_DIR_SOURCES}" PARENT_SCOPE) + + diff --git a/src/ross_files/Checkpoint.hh b/src/ross_files/Checkpoint.hh new file mode 100644 index 0000000..092ad0b --- /dev/null +++ b/src/ross_files/Checkpoint.hh @@ -0,0 +1,46 @@ +#pragma once + +#include +#include +#include + +// This could be done any number of ways faster/better than a list, +// but keep it for now. We can optimize this later, so long as we +// have the following interface. +class CheckpointQueue { + public: + void push_back(std::string newString) { queue.push_back(newString); } + void pop_back() { queue.pop_back(); } + std::string back() { return queue.back(); } + std::string front() { return queue.front(); } + void pop_front() { queue.pop_front(); } + + std::list queue; +}; + +class InputCheckpoint { + public: + std::string str() { return sss.str(); } + + template + InputCheckpoint& operator&(const TTT& ttt) { + sss.write(reinterpret_cast(&ttt), sizeof(TTT)); + return *this; + } + + std::stringstream sss; +}; + +class OutputCheckpoint { + public: + OutputCheckpoint(std::string newString) { + sss.str(newString); + } + template + OutputCheckpoint& operator&(TTT& ttt) { + sss.read(reinterpret_cast(&ttt), sizeof(TTT)); + return *this; + } + + std::stringstream sss; +}; diff --git a/src/ross_files/Funnel.hh b/src/ross_files/Funnel.hh new file mode 100644 index 0000000..c11a689 --- /dev/null +++ b/src/ross_files/Funnel.hh @@ -0,0 +1,560 @@ +#pragma once + +#include "actor.hh" +#include +#include + +#define MAX_TIME 100 //FIXME! +typedef std::size_t SpeciesTag; +typedef std::size_t ReactionTag; +typedef double Real; +typedef std::size_t SpeciesValue; +typedef Real (*RateFunction)(const std::vector& species); + +struct NoopMsg { +}; +struct RateMsg { + Real rate; +}; +struct ReactionMsg { + ReactionTag reaction; +}; +struct FunnelRateExchangeMsg { + ReactionTag reaction; + bool nrmMode; + Real rate; + Real normalContrib; +}; +struct SwitchModeMsg { + ReactionTag reaction; + bool nrmMode; +}; + +class NextReactionMethodCrucible : public LP { + public: + + ACTOR_DEF(); + + SLOT(firedMyself, NoopMsg); + void forward(firedMyself)(Time tnow, NoopMsg& msg, tw_bf*) { + ReactionMsg r; + r.reaction = _rxntag; + fire(0, r); + _countdown = newCountdown(); + fireNextReaction(tnow); + } + void commit(firedMyself)(Time tnow, NoopMsg& msg, tw_bf*) { + //std::cout << tnow << ": Fired Reaction " << _rxntag << "!\n"; + } + template + void checkpoint(firedMyself)(Checkpointer& ar, Time tnow, NoopMsg& msg, tw_bf*) { + ar & _countdown & _messageOutstanding & _lastUpdatedTime & _simpleRand; + } + + SLOT(updatedRate, RateMsg); + void forward(updatedRate)(Time tnow, RateMsg& msg, tw_bf*) { + if (_rate != 0) { + EVENT_CANCEL(_messageOutstanding); + _countdown -= countdownChange(tnow); + } + _rate = msg.rate; + if (_rate != 0) { + fireNextReaction(tnow); + } + } + template + void checkpoint(updatedRate)(Checkpointer& ar, Time tnow, RateMsg& msg, tw_bf*) { + ar & _messageOutstanding & _countdown & _rate & _lastUpdatedTime & _simpleRand; + } + + NextReactionMethodCrucible() { + _rate = 0; + _countdown = newCountdown(); + _messageOutstanding = PDES_NULL_EVENT; + } + + SIGNAL(fire,ReactionMsg); + + void setTag(const ReactionTag mytag) { + _rxntag = mytag; + } + + void setSeed(const int seed) { + _simpleRand.seed(seed); + } + + private: + ReactionTag _rxntag; + Real _countdown; + Real _rate; + std:: mt19937 _simpleRand; + PdesEvent _messageOutstanding; + Real _lastUpdatedTime; + // PdesLpid _self_lpid; + + void fireNextReaction(Time tnow) { + _lastUpdatedTime = tnow; + Time timeToFire = firingTime(_countdown,_rate); + NoopMsg msg; + _messageOutstanding = EVENT_SEND(twlp->gid, firedMyself, timeToFire, msg); + } + + Real newCountdown() { + //From Anderson, Algo 3 + //choose a uniform random number 0-1 + // std::random_device rd; //Will be used to obtain a seed for the random number engine + // std::mt19937 gen(rd()); //Standard mersenne_twister_engine seeded with rd() + std::uniform_real_distribution<> dis(0.0, 1.0); + Real uniform_random_number; + uniform_random_number = dis(_simpleRand); + + // std::cout << "Uniform: " << uniform_random_number << '\n'; + return log(1/uniform_random_number); + } + + Real countdownChange(Time tnow) { + return _rate*(tnow-_lastUpdatedTime); + } + + Real firingTime(Real countdown, Real rate) { + return countdown/rate; + } +}; + +extern +std::vector> stoic; + +class Funnel : public LP { + public: + + ACTOR_DEF(); + + Funnel() { + //Things needed to initiate the model + _rateNRM = 0; + _messageOutstanding = PDES_NULL_EVENT; + _numReactionsInBox = 0; + + + //tolerances the users can adjust + _tolerance = 0.003; + _sigma_max = 6; + _numReactionCutoff = 100; + _maxTimestep = MAX_TIME; + } + + + SLOT(firedReaction, ReactionMsg); + void forward(firedReaction)(Time tnow, ReactionMsg& msg, tw_bf* twbf) { + //update species + for (const auto& kv : stoic[msg.reaction]) { + auto iter = _indexFromSpecTag.find(kv.first); + if (iter != _indexFromSpecTag.end()) { + _species[iter->second] += kv.second; + } + } + updateSpeciesBasedOnRates(tnow); + bool allInBounds = checkRateBounds(tnow); + if (!allInBounds) { + interruptNextWakeup(tnow); + setupFutureEvents(tnow, true); + } + } + void backward(firedReaction)(Time tnow, ReactionMsg& msg, tw_bf*) { + } + void commit(firedReaction)(Time tnow, ReactionMsg& msg, tw_bf*) { + //std::cout << tnow << ": NRM " << msg.reaction << "->" << _rxntag << std::endl; + //printSpecies(tnow); + } + template + void checkpoint(firedReaction)(Checkpointer& ar, Time tnow, ReactionMsg& msg, tw_bf*) { + } + + + SLOT(changedRate, FunnelRateExchangeMsg); + void forward(changedRate)(Time tnow, FunnelRateExchangeMsg& msg, tw_bf* twbf) { + updateSpeciesBasedOnRates(tnow); + updateSavedRates(msg.reaction, msg.nrmMode, msg.rate, msg.normalContrib); + interruptNextWakeup(tnow); + //std::cout << tnow << ": %%% " << ((tnow-_lastRecomputeTime)>=_maxTimestep) << " " << checkRateBounds(tnow) << std::endl; + setupFutureEvents(tnow, (tnow-_lastRecomputeTime)>=_maxTimestep || !checkRateBounds(tnow)); + } + void backward(changedRate)(Time tnow, FunnelRateExchangeMsg& msg, tw_bf*) { + } + void commit(changedRate)(Time tnow, FunnelRateExchangeMsg& msg, tw_bf*) { + //std::cout << tnow << ": Changed rate for " << msg.reaction << "->" + // << _rxntag << ": " + // << (msg.nrmMode ? "NRM " : "SDE ") + // << msg.rate << " " + // << msg.normalContrib << std::endl; + //printSpecies(tnow); + } + template + void checkpoint(changedRate)(Checkpointer& ar, Time tnow, FunnelRateExchangeMsg& msg, tw_bf*) { + } + + SLOT(wakeup, NoopMsg); + void forward(wakeup)(Time tnow, NoopMsg& msg, tw_bf* twbf) { + // This gets called when we're sure that we're going to go + // out-of-bounds due to our tau condition, based on what we + // currently know about the SDE rates. + updateSpeciesBasedOnRates(tnow); + setupFutureEvents(tnow, (tnow-_lastRecomputeTime)>=_maxTimestep || !checkRateBounds(tnow)); + } + void commit(wakeup)(Time tnow, NoopMsg& msg, tw_bf*) { + //std::cout << tnow << ": Woke Reaction " << _rxntag << "!\n"; + //printSpecies(tnow); + } + template + void checkpoint(wakeup)(Checkpointer& ar, Time tnow, NoopMsg& msg, tw_bf*) { + } + + SLOT(changedMode, SwitchModeMsg); + void forward(changedMode)(Time tnow, SwitchModeMsg& msg, tw_bf*) { + updateSpeciesBasedOnRates(tnow); + int irxn = _indexFromRxnTag[msg.reaction]; + updateSavedRates(msg.reaction, msg.nrmMode, _rxnRate[irxn], _normalContrib[irxn]); + bool allInBounds = checkRateBounds(tnow); + if (!allInBounds) { + interruptNextWakeup(tnow); + setupFutureEvents(tnow, true); + } + } + void commit(changedMode)(Time tnow, SwitchModeMsg& msg, tw_bf*) { + //std::cout << tnow << ": Changed mode " << msg.reaction << "->" << _rxntag << ": " << msg.nrmMode << std::endl; + } + template + void checkpoint(changedMode)(Checkpointer& ar, Time tnow, SwitchModeMsg& msg, tw_bf*) { + } + + SLOT(printData, NoopMsg); + void forward(printData)(Time tnow, NoopMsg& msg, tw_bf*) { + } + void commit(printData)(Time tnow, NoopMsg& msg, tw_bf*) { + printSpecies(tnow); + } + + SIGNAL(updateRateNRM, RateMsg); + SIGNAL(updateRateFunnel, FunnelRateExchangeMsg); + SIGNAL(updateModeFunnel, SwitchModeMsg); + + void addSpecies(int tag, SpeciesValue initialValue, int trustRegion) { + int oldSize = _indexFromSpecTag.size(); + _indexFromSpecTag[tag] = oldSize; + _tagFromSpecIndex.push_back(tag); + _species.push_back(initialValue); + _lowerBound.push_back(0); + _upperBound.push_back(0); + _speciesRate.push_back(0); + _trustRegion.push_back(trustRegion); + } + + void setTag(int tag) { + _rxntag = tag; + } + + void setSeed(const int seed) { + _simpleRand.seed(seed); + } + + void setupSendToReaction(int otherRxnTag, tw_lpid destID) { + if (otherRxnTag != _rxntag) { + CONNECT(*this,Funnel::updateRateFunnel,destID,Funnel::changedRate); + CONNECT(*this,Funnel::updateModeFunnel,destID,Funnel::changedMode); + } + } + + void setupRecvFromReaction(int otherRxnTag) { + int index = _tagFromRxnIndex.size(); + _tagFromRxnIndex.push_back(otherRxnTag); + _indexFromRxnTag[otherRxnTag] = index; + _sqrtTimestep.push_back(0); + _normalContrib.push_back(0); + _rxnRate.push_back(0); + _rxnOverflow.push_back(0); + _rxnNrmMode.push_back(0); + if (otherRxnTag == _rxntag) { + _thisrxnindex = index; + } + } + + bool dependsOnSpecies(int specTag) { + return _indexFromSpecTag.find(specTag) != _indexFromSpecTag.end(); + } + + void beginReactions(Time tnow) { + _lastRecomputeTime = tnow; + _lastUpdatedTime = tnow; + + //recalculate rate + setupFutureEvents(tnow, true); + printSpecies(tnow); + } + + void printSpecies(Time tnow) { + std::cout << tnow << ": " + << "Reaction " << _rxntag << ": (wakeup " << _wakeupTime << ") (nrm=" << _rxnNrmMode[_thisrxnindex] << " rate=" << _rxnRate[_thisrxnindex] << ")" ; + for (int ii=0; ii<_tagFromSpecIndex.size(); ii++) { + std::cout << " {" << _tagFromSpecIndex[ii] << ": " + << "(" << _lowerBound[ii] << "<" + << _species[ii] + << "<" << _upperBound[ii] << ") " + << _speciesRate[ii] + << "}"; + } + std::cout << std::endl; + } + + public: + + RateFunction _rateFunction; //function to use for rate calculations + Real _tolerance; // tolerance for the tau independence condition + Real _sigma_max; // tolerance for switching between SDE and NRM + int _numReactionCutoff; // min number of reactions before we consider SDE + Real _maxTimestep; // maximum timestep allowed in SDE mode. + + + private: + std::unordered_map _indexFromSpecTag; // internal index from tag + std::vector _tagFromSpecIndex; // tag from internal index + std::vector _species; //species counts + std::vector _speciesRate; //first derivative of the species + std::vector _lowerBound; // lb of the tau independence region + std::vector _upperBound; // ub of the tau independence region + std::vector _trustRegion; //how far should we trust this species for a trust region? Think of this as a max for the ub/lb, independent of the derivative + + std::unordered_map _indexFromRxnTag; //internal index from tag + std::vector _tagFromRxnIndex; // tag from internal index + std::vector _sqrtTimestep; // sqrt of the timestep, used for SDE + std::vector _normalContrib; // sqrt(rate)*normal, used for SDE + std::vector _rxnOverflow; // Accumulator for SDE updates + std::vector _rxnRate; // rate for each reaction + std::vector _rxnNrmMode; // boolean, is the rxn in nrmMode + + Real _lastUpdatedTime; + Real _lastRecomputeTime; + Real _wakeupTime; + + ReactionTag _rxntag; + int _thisrxnindex; + + Real _rateNRM; + PdesEvent _messageOutstanding; + Real _numReactionsInBox; + + std:: mt19937 _simpleRand; + + Real recalculateRate(Time tnow) { + //recalculate rate + Real newRate = _rateFunction(_species); + _numReactionsInBox = _numReactionCutoff+1; + //compute gradient + std::vector gradient(_species.size()); + for (unsigned int ii=0; ii<_species.size(); ii++) { + std::vector modSpecies(_species); + modSpecies[ii]++; + Real gradRate = _rateFunction(modSpecies); + gradient[ii] = gradRate-newRate; + } + + int nonZeroGradCount = 0; + for (auto grad : gradient) { + if (grad!=0) { + nonZeroGradCount++; + } + } + //compute bounding box + for (unsigned int ii=0; ii<_species.size(); ii++) { + Real maxDelta = _trustRegion[ii]; + Real delta; + if (gradient[ii] == 0) { + delta = maxDelta; + } else { + delta = _tolerance*newRate/fabs(gradient[ii])/sqrt(nonZeroGradCount)/2; + if (delta > maxDelta) { + delta = maxDelta; + } + } + _lowerBound[ii] = _species[ii]-delta; + _upperBound[ii] = _species[ii]+delta; + + SpeciesTag specTag = _tagFromSpecIndex[ii]; + auto iter = stoic[_rxntag].find(specTag); + if (iter != stoic[_rxntag].end()) { + int sss = ((iter->second >= 0) ? iter->second : -iter->second); + if (sss) { + _numReactionsInBox = std::min(_numReactionsInBox, delta/sss); + } + } + } + _lastRecomputeTime = tnow; + return newRate; + } + + bool checkRateBounds(Time tnow) const { + //are we within our bounds + bool allWithinBounds=true; + for (unsigned int ii=0; ii<_species.size(); ii++) { + if (!(_lowerBound[ii] <= _species[ii] && _species[ii] <= _upperBound[ii])) { + allWithinBounds = false; + break; + } + } + return allWithinBounds; + } + + void updateSpeciesBasedOnRates(Time tnow) { + + // Use SDE euler method to update the rates for things in SDE + // mode. + if (tnow <= _lastUpdatedTime) { + return; + } + + for (unsigned int irxn=0; irxn<_tagFromRxnIndex.size(); irxn++) { + if (_rxnNrmMode[irxn]==false && _rxnRate[irxn] != 0) { + _rxnOverflow[irxn] += _rxnRate[irxn]*(tnow-_lastUpdatedTime); + if (_normalContrib[irxn] != 0) { + Real c_old = _sqrtTimestep[irxn]; + Real c_now = std::sqrt(tnow - _lastUpdatedTime + c_old * c_old); + _rxnOverflow[irxn] += _normalContrib[irxn]*(c_now-c_old); + _sqrtTimestep[irxn] = c_now; + } + if (_rxnOverflow[irxn] >= 1) { + int rxnCount = int(_rxnOverflow[irxn]); + for (const auto& kv : stoic[_tagFromRxnIndex[irxn]]) { + auto iter = _indexFromSpecTag.find(kv.first); + if (iter != _indexFromSpecTag.end()) { + _species[iter->second] += rxnCount*kv.second; + } + } + _rxnOverflow[irxn] -= rxnCount; + } + } + } + _lastUpdatedTime = tnow; + } + + void updateSavedRates(ReactionTag otherRxnTag, bool nrmMode, Real rate, Real normalContrib) { + int irxn = _indexFromRxnTag[otherRxnTag]; + _normalContrib[irxn] = normalContrib; + _sqrtTimestep[irxn] = 0; + Real delta_rate = rate-_rxnRate[irxn]; + _rxnRate[irxn] = rate; + _rxnNrmMode[irxn] = nrmMode; + if (delta_rate!=0) { + for (const auto& kv : stoic[otherRxnTag]) { + auto iter = _indexFromSpecTag.find(kv.first); + if (iter != _indexFromSpecTag.end()) { + _speciesRate[iter->second] += kv.second*delta_rate; + } + } + } + } + + void interruptNextWakeup(Time tnow) { + EVENT_CANCEL(_messageOutstanding); + } + + void scheduleNextWakeup(Time tnow, Real delta_time) { + NoopMsg msg; + _messageOutstanding = EVENT_SEND(twlp->gid, wakeup, delta_time, msg); + _wakeupTime = tnow+delta_time; + } + + void setupFutureEvents(Time tnow, bool recompute) { + + bool nrmMode; + Real delta_time; + Real nextRate; + if (recompute) { + nextRate = recalculateRate(tnow); + } else { + nextRate = _rxnRate[_thisrxnindex]; + } + getNextMode(nextRate, delta_time, nrmMode); + if (recompute) { + FunnelRateExchangeMsg newRateMsg; + newRateMsg.reaction = _rxntag; + newRateMsg.rate = nextRate; + newRateMsg.nrmMode = nrmMode; + + std::normal_distribution<> dis(0,1); + Real unit_normal_random = dis(_simpleRand); + newRateMsg.normalContrib = unit_normal_random * std::sqrt(nextRate); + + updateSavedRates(newRateMsg.reaction, newRateMsg.nrmMode, newRateMsg.rate, newRateMsg.normalContrib); + updateRateFunnel(0, newRateMsg); + } + Real newRateNRM = (nrmMode==true ? nextRate : 0); + if (newRateNRM != _rateNRM) { + RateMsg newRateMsg; + newRateMsg.rate = newRateNRM; + updateRateNRM(0, newRateMsg); + _rateNRM = newRateNRM; + + if (!recompute) { + SwitchModeMsg newSwitchMsg; + newSwitchMsg.reaction = _rxntag; + newSwitchMsg.nrmMode = nrmMode; + updateSavedRates(_rxntag, nrmMode, nextRate, _normalContrib[_thisrxnindex]); + updateModeFunnel(0, newSwitchMsg); + } + } + scheduleNextWakeup(tnow, delta_time); + } + + void getNextMode(Real nextRate, Real& delta_time, bool& nextNrmMode) const { + + delta_time = findWakeupInterval(nextRate); + if (_numReactionsInBox > _numReactionCutoff && nextRate*(delta_time+_sqrtTimestep[_thisrxnindex]*_sqrtTimestep[_thisrxnindex]) > _sigma_max*_sigma_max) { + nextNrmMode = false; + } else { + nextNrmMode = true; + } + } + + Real findWakeupInterval(Real nextRate) const { + // find when species will cross boundary + Real minTime = _maxTimestep; + Real delta_rate = nextRate-(_rxnRate[_thisrxnindex]); + + for (unsigned int ispec=0; ispec<_species.size(); ispec++) { + Real thisRate = _speciesRate[ispec]; + if (delta_rate != 0) { + auto iter=stoic[_rxntag].find(_tagFromSpecIndex[ispec]); + if (iter != stoic[_rxntag].end()) { + thisRate += delta_rate*iter->second; + } + } + if (thisRate > 0) { + minTime = std::min(minTime,(_upperBound[ispec]+1 - _species[ispec]) / thisRate); + } else if (thisRate < 0) { + minTime = std::min(minTime,(_lowerBound[ispec]-1 - _species[ispec]) / thisRate); + } else { + //do nothing + } + //std::cout << "%%% " << thisRate << " " << minTime << std::endl; + } + return minTime; + } +}; + + +class Heartbeat : public LP { + public: + + ACTOR_DEF(); + + void setInterval(double interval) { _interval = interval; } + + SLOT(activated, NoopMsg); + void forward(activated)(Time tnow, NoopMsg& msg, tw_bf*) { + activate(_interval, msg); + } + + SIGNAL(activate, NoopMsg); + + private: + double _interval; +}; diff --git a/src/ross_files/actor.cc b/src/ross_files/actor.cc new file mode 100644 index 0000000..1d97998 --- /dev/null +++ b/src/ross_files/actor.cc @@ -0,0 +1,13 @@ + +#include "actor.hh" +#include + +void my_event_cancel(PdesEvent evtptr) { + if (evtptr) { + tw_event_rescind(evtptr); + } +} + +tw_lpid BaseLP::getID() const { + return twlp->gid; +} diff --git a/src/ross_files/actor.hh b/src/ross_files/actor.hh new file mode 100644 index 0000000..47978ab --- /dev/null +++ b/src/ross_files/actor.hh @@ -0,0 +1,138 @@ +#pragma once + +#if defined(WCS_HAS_ROSS) +#include +#include +#include +#include +// #include "wcs-ross-bf.hpp" +#include "simtime.hh" +#include "Checkpoint.hh" + +//// Classes that make the whole API work + +struct EventData; //will be defined by a perl script. +typedef std::size_t Tag; +typedef tw_lpid PdesLpid; +typedef tw_event* PdesEvent; +#define PDES_NULL_EVENT NULL + +class BaseLP { + public: + virtual ~BaseLP() {} + + struct ObserverData { + tw_lpid dest; + Tag slot; + }; + tw_lp* twlp; //this has to be initialized properly! TODO + + tw_lpid getID() const; +}; + +extern CheckpointQueue checkpointQueue; + +template +class LP : public BaseLP { + public: + typedef void (Dependent::*MethodPtr)(Time tnow, void*, tw_bf*); + static void init_lp(void *state, tw_lp *me); + static void forward_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp); + static void backward_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp); + static void commit_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp); + static void final_lp(void *state, tw_lp *me); + static MethodPtr forwardDispatch[]; + static MethodPtr backwardDispatch[]; + static MethodPtr commitDispatch[]; + + template + void forwardWithCheckpoint(Time tnow, typename SlotType::MsgType& msg, tw_bf* twbf) { + InputCheckpoint ic; + static_cast(this)->Dependent::template checkpointFull(ic, tnow, msg, twbf); + static_cast(this)->Dependent::template forwardFull(tnow, msg, twbf); + checkpointQueue.push_back(ic.str()); + } + template + void backwardWithCheckpoint(Time tnow, typename SlotType::MsgType& msg, tw_bf* twbf) { + static_cast(this)->Dependent::template backwardFull(tnow, msg, twbf); + OutputCheckpoint oc(checkpointQueue.back()); + static_cast(this)->Dependent::template checkpointFull(oc, tnow, msg, twbf); + checkpointQueue.pop_back(); + } + template + void commitWithCheckpoint(Time tnow, typename SlotType::MsgType& msg, tw_bf* twbf) { + static_cast(this)->Dependent::template commitFull(tnow, msg, twbf); + checkpointQueue.pop_front(); + } +}; + +#define ACTOR_DEF() \ + template \ + void forwardFull(Time tnow, typename SlotType::MsgType& msg, tw_bf*); \ + template \ + void backwardFull(Time tnow, typename SlotType::MsgType& msg, tw_bf*) {} \ + template \ + void commitFull(Time tnow, typename SlotType::MsgType& msg, tw_bf*) {} \ + template \ + inline void checkpointFull(Archiver& ar, Time tnow, typename SlotType::MsgType& msg, tw_bf*); + + +#define TAGIFY(SlotType) SlotType::tag + +#define forward(SlotType) forward_##SlotType +#define backward(SlotType) backward_##SlotType +#define commit(SlotType) commit_##SlotType +#define checkpoint(SlotType) checkpoint_##SlotType + +#define SLOT(name, msgType) \ +class name { \ + public: \ + typedef msgType MsgType; \ + static Tag tag; \ +} + +#define SIGNAL(name, type) \ + class name##_signal_class { public: typedef type MsgType; }; \ + void name(Time trecv, type& msg) { \ + /*for all LPs attached to this signal*/ \ + for (auto&& observerData : name##_listeners) { \ + my_event_send(observerData.dest, trecv, twlp, observerData.slot, msg); \ + } \ + } \ + std::list name##_listeners + +#define CONNECT(sendingLP, signal, dest, slot) signal_add_listener((sendingLP).signal##_listeners,dest) + +template +inline +typename std::enable_if::value>::type +signal_add_listener(std::list& listeners, const BaseLP& dest) { + tw_lpid destID = dest.getID(); + signal_add_listener(listeners, destID); +} + +template +inline +typename std::enable_if::value>::type +signal_add_listener(std::list& listeners, const tw_lpid destID) { + listeners.push_back({destID,TAGIFY(SlotType)}); +} + +#define EVENT_SEND(lpid, slotClass, msgTime, msg) my_event_send(lpid, msgTime, twlp, msg) + +#define EVENT_CANCEL(event) my_event_cancel(event) + +template +PdesEvent my_event_send(tw_lpid dest, Time trecv, tw_lp* twlp, typename SlotType::MsgType& msg) +{ + return my_event_send(dest, trecv, twlp, TAGIFY(SlotType), msg); +} + +template +PdesEvent my_event_send(tw_lpid dest, Time trecv, tw_lp* twlp, Tag tag, TTT& msg); + +void my_event_cancel(PdesEvent evtptr); + +std::size_t messageSize(); + +#endif // defined(WCS_HAS_ROSS) \ No newline at end of file diff --git a/src/ross_files/generateInfrastructure.py b/src/ross_files/generateInfrastructure.py new file mode 100644 index 0000000..6a968f4 --- /dev/null +++ b/src/ross_files/generateInfrastructure.py @@ -0,0 +1,188 @@ +#!/usr/bin/env python + +import re +from collections import OrderedDict + +class Slot: + def __init__(self, name, msg): + self.name = name + self.msg = msg + self.modes = set() + def hasMode(self, mode): + return mode in self.modes + def addMode(self, mode): + self.modes.add(mode) + +class Actor: + def __init__(self, name): + self.name = name + self.slots = {} + def addSlot(self,slot): + self.slots[slot.name] = slot + def getSlot(self,slotName): + return self.slots[slotName] + def getSlots(self): + return self.slots.values() + +def main(): + import sys + + filenames = sys.argv[1:] + + actors = OrderedDict() + messages = set() + for filename in filenames: + for line in open(filename, "r"): + mmm = re.search(r'public LP<([A-Za-z0-9_]+)>',line) + if mmm: + actorName = mmm.group(1) + actors[actorName] = Actor(actorName) + mmm = re.search(r'SLOT\(\s*([A-Za-z0-9_]+)\s*,\s*([A-Za-z0-9_]+)\s*\)',line) + if mmm: + slotName = mmm.group(1) + msgName = mmm.group(2) + actors[actorName].addSlot(Slot(slotName,msgName)) + for eventType in ['forward','backward','commit','checkpoint']: + mmm = re.search("void "+eventType+r'\(([A-Za-z0-9_]+)\)',line) + if mmm: + slotName = mmm.group(1) + actors[actorName].getSlot(slotName).addMode(eventType) + + print('#include "actor.hh"') + print('#include "Checkpoint.hh"') + + for filename in filenames: + print('#include "%s"' % filename) + + print(''' +struct EventData { + Tag tag; + union { + char msg;''') + allMessages = set() + for actor in actors.values(): + for slot in actor.getSlots(): + allMessages.add(slot.msg) + ii=0 + for msg in allMessages: + print(" %s msg%d;" % (msg, ii)) + ii += 1 + print(''' + }; +}; + +std::size_t messageSize() { + return sizeof(EventData); +} + +CheckpointQueue checkpointQueue; + +template +PdesEvent my_event_send(tw_lpid dest, Time trecv, tw_lp* twlp, Tag tag, TTT& msg) +{ + tw_stime recv_ts = TW_STIME_ADD(tw_now(twlp), trecv); + if (recv_ts >= g_tw_ts_end) { + return NULL; + } + tw_event* evtptr = tw_event_new(dest, trecv, twlp); + EventData* data = static_cast(tw_event_data(evtptr)); + data->tag = tag; + memcpy(&data->msg, &msg, sizeof(msg)); + tw_event_send(evtptr); + return evtptr; +}''') + for msg in allMessages: + print("template PdesEvent my_event_send(tw_lpid dest, Time trecv, tw_lp* twlp, Tag tag, %s& msg);" % msg) + for actor in actors.values(): + for slot in actor.getSlots(): + ttt = { + "actor" : actor.name, + "slot" : slot.name, + "msg" : slot.msg, + } + for eventType in "forward","backward","commit": + ttt['event'] = eventType + if slot.hasMode(eventType): + print("""template <> +inline void %(actor)s::%(event)sFull<%(actor)s::%(slot)s>(Time tnow, %(msg)s& msg, tw_bf* twbf) { %(event)s(%(slot)s)(tnow,msg,twbf); }""" % ttt) + if slot.hasMode('checkpoint'): + print("""template <> +inline void %(actor)s::checkpointFull<%(actor)s::%(slot)s>(InputCheckpoint& ar, Time tnow, %(msg)s& msg, tw_bf* twbf) { checkpoint(%(slot)s)(ar,tnow,msg,twbf); } +template <> +inline void %(actor)s::checkpointFull<%(actor)s::%(slot)s>(OutputCheckpoint& ar, Time tnow, %(msg)s& msg, tw_bf* twbf) { checkpoint(%(slot)s)(ar,tnow,msg,twbf); } +""" % ttt) + for actor in actors.values(): + for (ii,slot) in enumerate(actor.getSlots()): + print("Tag %s::%s::tag = %d;" % (actor.name,slot.name,ii)) + + for actor in actors.values(): + for eventType in ["forward","backward","commit"]: + print("template<> %s::MethodPtr LP<%s>::%sDispatch[] = {" % (actor.name, actor.name,eventType)) + for slot in actor.getSlots(): + ttt = { + "actor" : actor.name, + "slot" : slot.name, + "event" : eventType, + } + if slot.hasMode('checkpoint'): + method = '%(actor)s::%(event)sWithCheckpoint<%(actor)s::%(slot)s>' % ttt + else: + method = '%(actor)s::%(event)sFull<%(actor)s::%(slot)s>' % ttt + ttt['method'] = method + print(" reinterpret_cast<%(actor)s::MethodPtr>(&%(method)s)," % ttt) + print("0};") + + print(""" +///This stuff needs EventData to be defined, but doesn't need specific generation + +template +void LP::init_lp(void *state, tw_lp *thislp) { + Dependent* derivedObj = new (reinterpret_cast(state)) Dependent(); + BaseLP *lpptr = static_cast(derivedObj); + lpptr->twlp = thislp; +} + +template +void LP::forward_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp) { + Time tnow = tw_now(thislp); + EventData *e = reinterpret_cast(evt); + Dependent* derivedObj = reinterpret_cast(state); + (derivedObj->*forwardDispatch[e->tag])(tnow, &e->msg, bf); +} + +template +void LP::backward_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp) { + Time tnow = tw_now(thislp); + EventData *e = reinterpret_cast(evt); + Dependent* derivedObj = reinterpret_cast(state); + (derivedObj->*backwardDispatch[e->tag])(tnow, &e->msg, bf); +} + +template +void LP::commit_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp) { + Time tnow = tw_now(thislp); + EventData *e = reinterpret_cast(evt); + Dependent* derivedObj = reinterpret_cast(state); + (derivedObj->*commitDispatch[e->tag])(tnow, &e->msg, bf); +} + +template +void LP::final_lp(void *state, tw_lp *) { + Dependent* derivedObj = reinterpret_cast(state); + derivedObj->~Dependent(); +} + +""") + for actor in actors.values(): + print(""" +template void LP<%(actor)s>::forward_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp); +template void LP<%(actor)s>::backward_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp); +template void LP<%(actor)s>::commit_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp); +template void LP<%(actor)s>::final_lp(void *state, tw_lp *me); +template void LP<%(actor)s>::init_lp(void *state, tw_lp *me); + +""" % {"actor" : actor.name}) + + +if __name__=="__main__": + main() diff --git a/src/ross_files/simtime.hh b/src/ross_files/simtime.hh new file mode 100644 index 0000000..6a5d415 --- /dev/null +++ b/src/ross_files/simtime.hh @@ -0,0 +1,10 @@ +#pragma once + +#ifndef KMCSIM_TYPES_ONLY +#include +#include +// #include "wcs-ross-bf.hpp" +#endif + +typedef tw_stime Time; + From ebd6d1adac47bde35bdf135917fb8d7acf638f17 Mon Sep 17 00:00:00 2001 From: Konstantia Georgouli Date: Fri, 27 Aug 2021 08:51:55 -0700 Subject: [PATCH 04/14] Generate actor.gen.cc --- CMakeLists.txt | 2 + src/ross_files/CMakeLists.txt | 51 ++++++++++++++++++------ src/ross_files/actor.hh | 3 -- src/ross_files/generateInfrastructure.py | 0 4 files changed, 41 insertions(+), 15 deletions(-) mode change 100644 => 100755 src/ross_files/generateInfrastructure.py diff --git a/CMakeLists.txt b/CMakeLists.txt index 899f545..ab1cd41 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -430,6 +430,8 @@ if (WCS_HAS_PROTOBUF) endif (WCS_HAS_OPENMP) endif (WCS_HAS_PROTOBUF) +add_dependencies(wcs wcs_actorgen) + target_include_directories(wcs PUBLIC $ $ diff --git a/src/ross_files/CMakeLists.txt b/src/ross_files/CMakeLists.txt index e8c70ed..8c1fe2b 100644 --- a/src/ross_files/CMakeLists.txt +++ b/src/ross_files/CMakeLists.txt @@ -5,7 +5,7 @@ # ) # find_package(Python COMPONENTS Interpreter) -ADD_CUSTOM_COMMAND(OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/actor.gen.cc +add_custom_command(OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/actor.gen.cc COMMAND ${Python_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/generateInfrastructure.py ${actor_srcs} > ${CMAKE_CURRENT_BINARY_DIR}/actor.gen.cc MAIN_DEPENDENCY generateInfrastructure.py DEPENDS generateInfrastructure.py ${actor_srcs} @@ -14,23 +14,50 @@ SET_PROPERTY(SOURCE ${CMAKE_CURRENT_BINARY_DIR}/actor.gen.cc PROPERTY INCLUDE_DIRECTORIES ${CMAKE_CURRENT_SOURCE_DIR}) +# SET(my_srcs +# ${actor_srcs} +# actor.hh +# actor.cc +# simtime.hh +# Checkpoint.hh +# ) +# # copy .hh and .cc fiels in the folder with actor.gen.cc +# foreach( file_i ${my_srcs}) +# add_custom_command( +# OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/${file_i} +# COMMAND ${CMAKE_COMMAND} -E copy ${CMAKE_CURRENT_SOURCE_DIR}/${file_i} ${CMAKE_CURRENT_BINARY_DIR}/${file_i} +# ) +# endforeach( file_i ) + +add_custom_target(wcs_actorgen + DEPENDS ${CMAKE_CURRENT_BINARY_DIR}/actor.gen.cc) + # DEPENDS ${CMAKE_CURRENT_BINARY_DIR}/actor.gen.cc ${CMAKE_CURRENT_BINARY_DIR}/${my_srcs}) + +# add_dependencies(wcs wcs_actorgen) # Add the header and source files for this directory -# set_full_path(THIS_DIR_HEADERS -# actor.hh -# Checkpoint.hh -# Funnel.hh -# simtime.hh -# ) +set_full_path(THIS_DIR_HEADERS + actor.hh + Checkpoint.hh + Funnel.hh + simtime.hh + ) -# set_full_path(THIS_DIR_SOURCES -# actor.cc -# ) +set_full_path(THIS_DIR_SOURCES + actor.cc + ) # # Add the subdirectories +set_source_files_properties( + ${CMAKE_CURRENT_BINARY_DIR}/actor.gen.cc + PROPERTIES GENERATED TRUE +) + +# set(THIS_DIR_SOURCES "${THIS_DIR_SOURCES}" "${CMAKE_CURRENT_BINARY_DIR}/actor.gen.cc") + # # Propagate the files up the tree -# set(WCS_HEADERS "${WCS_HEADERS}" "${THIS_DIR_HEADERS}" PARENT_SCOPE) -# set(WCS_SOURCES "${WCS_SOURCES}" "${THIS_DIR_SOURCES}" PARENT_SCOPE) +set(WCS_HEADERS "${WCS_HEADERS}" "${THIS_DIR_HEADERS}" PARENT_SCOPE) +set(WCS_SOURCES "${WCS_SOURCES}" "${THIS_DIR_SOURCES}" PARENT_SCOPE) diff --git a/src/ross_files/actor.hh b/src/ross_files/actor.hh index 47978ab..5dadc26 100644 --- a/src/ross_files/actor.hh +++ b/src/ross_files/actor.hh @@ -1,6 +1,5 @@ #pragma once -#if defined(WCS_HAS_ROSS) #include #include #include @@ -134,5 +133,3 @@ PdesEvent my_event_send(tw_lpid dest, Time trecv, tw_lp* twlp, Tag tag, TTT& msg void my_event_cancel(PdesEvent evtptr); std::size_t messageSize(); - -#endif // defined(WCS_HAS_ROSS) \ No newline at end of file diff --git a/src/ross_files/generateInfrastructure.py b/src/ross_files/generateInfrastructure.py old mode 100644 new mode 100755 From b095b545907578ce9d916bdbdb3811c6c9e8f89a Mon Sep 17 00:00:00 2001 From: Konstantia Georgouli Date: Fri, 27 Aug 2021 13:56:14 -0700 Subject: [PATCH 05/14] Fix linking with ROSS --- CMakeLists.txt | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index ab1cd41..dd50fde 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -457,6 +457,10 @@ if (WCS_HAS_SBML) # target_include_directories(wcs SYSTEM PUBLIC ${SBML_INCLUDE_DIRS}) endif (WCS_HAS_SBML) +if (WCS_HAS_ROSS) + target_link_libraries(wcs PUBLIC ROSS::ROSS) +endif (WCS_HAS_ROSS) + if (WCS_HAS_EXPRTK) add_dependencies(wcs EXPRTK-download) target_include_directories(wcs PUBLIC @@ -630,7 +634,7 @@ if (WCS_HAS_ROSS) $ $ $) - target_link_libraries(hybrid-bin PRIVATE ROSS::ROSS MPI::MPI_CXX wcs ${LIB_FILESYSTEM}) + target_link_libraries(hybrid-bin PRIVATE MPI::MPI_CXX wcs ${LIB_FILESYSTEM}) set_target_properties(hybrid-bin PROPERTIES OUTPUT_NAME hybrid) list(APPEND WCS_EXEC_TARGETS hybrid-bin) endif (WCS_HAS_ROSS) From 647343bc588b030933782839332f40af0f8ef97a Mon Sep 17 00:00:00 2001 From: Robert Blake Date: Wed, 1 Sep 2021 16:33:44 -0700 Subject: [PATCH 06/14] Code builds now. --- CMakeLists.txt | 18 ++-------- src/CMakeLists.txt | 1 - src/ross_files/CMakeLists.txt | 59 +++++++++++++++------------------ src/ross_files/Funnel.hh | 6 ++-- src/{ => ross_files}/hybrid.cpp | 0 src/{ => ross_files}/hybrid.hpp | 0 6 files changed, 32 insertions(+), 52 deletions(-) rename src/{ => ross_files}/hybrid.cpp (100%) rename src/{ => ross_files}/hybrid.hpp (100%) diff --git a/CMakeLists.txt b/CMakeLists.txt index dd50fde..e1f519e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -430,8 +430,6 @@ if (WCS_HAS_PROTOBUF) endif (WCS_HAS_OPENMP) endif (WCS_HAS_PROTOBUF) -add_dependencies(wcs wcs_actorgen) - target_include_directories(wcs PUBLIC $ $ @@ -457,10 +455,6 @@ if (WCS_HAS_SBML) # target_include_directories(wcs SYSTEM PUBLIC ${SBML_INCLUDE_DIRS}) endif (WCS_HAS_SBML) -if (WCS_HAS_ROSS) - target_link_libraries(wcs PUBLIC ROSS::ROSS) -endif (WCS_HAS_ROSS) - if (WCS_HAS_EXPRTK) add_dependencies(wcs EXPRTK-download) target_include_directories(wcs PUBLIC @@ -621,21 +615,13 @@ if (WCS_HAS_ROSS) $ $ $) - target_link_libraries(des-bin PRIVATE ROSS::ROSS MPI::MPI_CXX wcs ${LIB_FILESYSTEM}) + target_link_libraries(des-bin PRIVATE wcs ROSS::ROSS MPI::MPI_CXX ${LIB_FILESYSTEM}) set_target_properties(des-bin PROPERTIES OUTPUT_NAME des) list(APPEND WCS_EXEC_TARGETS des-bin) endif (WCS_HAS_ROSS) # add executable hybrid if (WCS_HAS_ROSS) - add_executable( hybrid-bin src/hybrid.cpp ) - target_include_directories(hybrid-bin SYSTEM PUBLIC ${ROSS_INCLUDE_DIRS}) - target_include_directories(hybrid-bin PUBLIC - $ - $ - $) - target_link_libraries(hybrid-bin PRIVATE MPI::MPI_CXX wcs ${LIB_FILESYSTEM}) - set_target_properties(hybrid-bin PROPERTIES OUTPUT_NAME hybrid) list(APPEND WCS_EXEC_TARGETS hybrid-bin) endif (WCS_HAS_ROSS) @@ -816,7 +802,7 @@ install( "${PROJECT_SOURCE_DIR}/src/wcs_types.hpp" "${PROJECT_SOURCE_DIR}/src/bgl.hpp" "${PROJECT_SOURCE_DIR}/src/des.hpp" - "${PROJECT_SOURCE_DIR}/src/hybrid.hpp" + "${PROJECT_SOURCE_DIR}/src/ross_files/hybrid.hpp" "${PROJECT_SOURCE_DIR}/src/wcs-ross-bf.hpp" DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}") diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ab109fc..da4a958 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -3,7 +3,6 @@ set_full_path(THIS_DIR_HEADERS bgl.hpp wcs_types.hpp des.hpp - hybrid.hpp wcs-ross-bf.hpp ) diff --git a/src/ross_files/CMakeLists.txt b/src/ross_files/CMakeLists.txt index 8c1fe2b..1c83b3c 100644 --- a/src/ross_files/CMakeLists.txt +++ b/src/ross_files/CMakeLists.txt @@ -5,54 +5,49 @@ # ) # find_package(Python COMPONENTS Interpreter) +set(actor_srcs ${CMAKE_CURRENT_SOURCE_DIR}/Funnel.hh) + add_custom_command(OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/actor.gen.cc COMMAND ${Python_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/generateInfrastructure.py ${actor_srcs} > ${CMAKE_CURRENT_BINARY_DIR}/actor.gen.cc MAIN_DEPENDENCY generateInfrastructure.py DEPENDS generateInfrastructure.py ${actor_srcs} VERBATIM) -SET_PROPERTY(SOURCE ${CMAKE_CURRENT_BINARY_DIR}/actor.gen.cc - PROPERTY INCLUDE_DIRECTORIES ${CMAKE_CURRENT_SOURCE_DIR}) +set_source_files_properties( + ${CMAKE_CURRENT_BINARY_DIR}/actor.gen.cc + PROPERTIES GENERATED TRUE +) -# SET(my_srcs -# ${actor_srcs} -# actor.hh -# actor.cc -# simtime.hh -# Checkpoint.hh -# ) +if (WCS_HAS_ROSS) + add_executable(hybrid-bin + hybrid.cpp + hybrid.hpp + actor.hh + Checkpoint.hh + simtime.hh + actor.cc + ${actor_srcs} + ${CMAKE_CURRENT_BINARY_DIR}/actor.gen.cc + ) + target_link_libraries(hybrid-bin PUBLIC ROSS::ROSS MPI::MPI_CXX wcs ${LIB_FILESYSTEM}) + target_include_directories(hybrid-bin SYSTEM PUBLIC ${ROSS_INCLUDE_DIRS}) + target_include_directories(hybrid-bin PUBLIC + $ + $ + $ + $) + + set_target_properties(hybrid-bin PROPERTIES OUTPUT_NAME hybrid) +endif() -# # copy .hh and .cc fiels in the folder with actor.gen.cc -# foreach( file_i ${my_srcs}) -# add_custom_command( -# OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/${file_i} -# COMMAND ${CMAKE_COMMAND} -E copy ${CMAKE_CURRENT_SOURCE_DIR}/${file_i} ${CMAKE_CURRENT_BINARY_DIR}/${file_i} -# ) -# endforeach( file_i ) - -add_custom_target(wcs_actorgen - DEPENDS ${CMAKE_CURRENT_BINARY_DIR}/actor.gen.cc) - # DEPENDS ${CMAKE_CURRENT_BINARY_DIR}/actor.gen.cc ${CMAKE_CURRENT_BINARY_DIR}/${my_srcs}) - -# add_dependencies(wcs wcs_actorgen) -# Add the header and source files for this directory set_full_path(THIS_DIR_HEADERS - actor.hh - Checkpoint.hh - Funnel.hh - simtime.hh ) set_full_path(THIS_DIR_SOURCES - actor.cc ) # # Add the subdirectories -set_source_files_properties( - ${CMAKE_CURRENT_BINARY_DIR}/actor.gen.cc - PROPERTIES GENERATED TRUE -) # set(THIS_DIR_SOURCES "${THIS_DIR_SOURCES}" "${CMAKE_CURRENT_BINARY_DIR}/actor.gen.cc") diff --git a/src/ross_files/Funnel.hh b/src/ross_files/Funnel.hh index c11a689..3283f7b 100644 --- a/src/ross_files/Funnel.hh +++ b/src/ross_files/Funnel.hh @@ -257,14 +257,14 @@ class Funnel : public LP { _simpleRand.seed(seed); } - void setupSendToReaction(int otherRxnTag, tw_lpid destID) { + void setupSendToReaction(ReactionTag otherRxnTag, tw_lpid destID) { if (otherRxnTag != _rxntag) { CONNECT(*this,Funnel::updateRateFunnel,destID,Funnel::changedRate); CONNECT(*this,Funnel::updateModeFunnel,destID,Funnel::changedMode); } } - void setupRecvFromReaction(int otherRxnTag) { + void setupRecvFromReaction(ReactionTag otherRxnTag) { int index = _tagFromRxnIndex.size(); _tagFromRxnIndex.push_back(otherRxnTag); _indexFromRxnTag[otherRxnTag] = index; @@ -294,7 +294,7 @@ class Funnel : public LP { void printSpecies(Time tnow) { std::cout << tnow << ": " << "Reaction " << _rxntag << ": (wakeup " << _wakeupTime << ") (nrm=" << _rxnNrmMode[_thisrxnindex] << " rate=" << _rxnRate[_thisrxnindex] << ")" ; - for (int ii=0; ii<_tagFromSpecIndex.size(); ii++) { + for (std::size_t ii=0; ii<_tagFromSpecIndex.size(); ii++) { std::cout << " {" << _tagFromSpecIndex[ii] << ": " << "(" << _lowerBound[ii] << "<" << _species[ii] diff --git a/src/hybrid.cpp b/src/ross_files/hybrid.cpp similarity index 100% rename from src/hybrid.cpp rename to src/ross_files/hybrid.cpp diff --git a/src/hybrid.hpp b/src/ross_files/hybrid.hpp similarity index 100% rename from src/hybrid.hpp rename to src/ross_files/hybrid.hpp From 71279341fc2518c66220b9a8628e131d313b1bdd Mon Sep 17 00:00:00 2001 From: Konstantia Georgouli Date: Thu, 2 Sep 2021 07:32:37 -0700 Subject: [PATCH 07/14] Rename ross_files folder to hybrid --- CMakeLists.txt | 6 +++--- src/CMakeLists.txt | 2 +- src/{ross_files => hybrid}/CMakeLists.txt | 2 +- src/{ross_files => hybrid}/Checkpoint.hh | 0 src/{ross_files => hybrid}/Funnel.hh | 0 src/{ross_files => hybrid}/actor.cc | 0 src/{ross_files => hybrid}/actor.hh | 0 src/{ross_files => hybrid}/generateInfrastructure.py | 0 src/{ross_files => hybrid}/hybrid.cpp | 4 ++-- src/{ross_files => hybrid}/hybrid.hpp | 0 src/{ross_files => hybrid}/simtime.hh | 0 11 files changed, 7 insertions(+), 7 deletions(-) rename src/{ross_files => hybrid}/CMakeLists.txt (98%) rename src/{ross_files => hybrid}/Checkpoint.hh (100%) rename src/{ross_files => hybrid}/Funnel.hh (100%) rename src/{ross_files => hybrid}/actor.cc (100%) rename src/{ross_files => hybrid}/actor.hh (100%) rename src/{ross_files => hybrid}/generateInfrastructure.py (100%) rename src/{ross_files => hybrid}/hybrid.cpp (99%) rename src/{ross_files => hybrid}/hybrid.hpp (100%) rename src/{ross_files => hybrid}/simtime.hh (100%) diff --git a/CMakeLists.txt b/CMakeLists.txt index e1f519e..61bb9c3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -718,7 +718,7 @@ set(INCLUDE_INSTALL_DIRS "${CMAKE_SOURCE_DIR}/src/proto" "${CMAKE_SOURCE_DIR}/src/reaction_network" "${CMAKE_SOURCE_DIR}/src/sim_methods" - "${CMAKE_SOURCE_DIR}/src/ross_files" + "${CMAKE_SOURCE_DIR}/src/hybrid" "${CMAKE_SOURCE_DIR}/src/utils") set(LIB_INSTALL_DIR "${CMAKE_BINARY_DIR}") set(EXTRA_CMAKE_MODULE_DIR "${CMAKE_SOURCE_DIR}/cmake/modules") @@ -794,7 +794,7 @@ install( DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}" FILES_MATCHING PATTERN "*.h" PATTERN "*.hpp") install( - DIRECTORY "${PROJECT_SOURCE_DIR}/src/ross_files" + DIRECTORY "${PROJECT_SOURCE_DIR}/src/hybrid" DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}" FILES_MATCHING PATTERN "*.h" PATTERN "*.hpp" PATTERN "*.hh") install( @@ -802,7 +802,7 @@ install( "${PROJECT_SOURCE_DIR}/src/wcs_types.hpp" "${PROJECT_SOURCE_DIR}/src/bgl.hpp" "${PROJECT_SOURCE_DIR}/src/des.hpp" - "${PROJECT_SOURCE_DIR}/src/ross_files/hybrid.hpp" + "${PROJECT_SOURCE_DIR}/src/hybrid/hybrid.hpp" "${PROJECT_SOURCE_DIR}/src/wcs-ross-bf.hpp" DESTINATION "${CMAKE_INSTALL_INCLUDEDIR}") diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index da4a958..e55c471 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -16,7 +16,7 @@ add_subdirectory(proto) add_subdirectory(reaction_network) add_subdirectory(sim_methods) add_subdirectory(utils) -add_subdirectory(ross_files) +add_subdirectory(hybrid) # Propagate the files up the tree set(WCS_HEADERS "${WCS_HEADERS}" "${THIS_DIR_HEADERS}" PARENT_SCOPE) diff --git a/src/ross_files/CMakeLists.txt b/src/hybrid/CMakeLists.txt similarity index 98% rename from src/ross_files/CMakeLists.txt rename to src/hybrid/CMakeLists.txt index 1c83b3c..51c3837 100644 --- a/src/ross_files/CMakeLists.txt +++ b/src/hybrid/CMakeLists.txt @@ -1,5 +1,5 @@ # blt_add_executable( -# NAME ross_files +# NAME hybrid # SOURCES actor.cc # DEPENDS_ON ROSS mpi # ) diff --git a/src/ross_files/Checkpoint.hh b/src/hybrid/Checkpoint.hh similarity index 100% rename from src/ross_files/Checkpoint.hh rename to src/hybrid/Checkpoint.hh diff --git a/src/ross_files/Funnel.hh b/src/hybrid/Funnel.hh similarity index 100% rename from src/ross_files/Funnel.hh rename to src/hybrid/Funnel.hh diff --git a/src/ross_files/actor.cc b/src/hybrid/actor.cc similarity index 100% rename from src/ross_files/actor.cc rename to src/hybrid/actor.cc diff --git a/src/ross_files/actor.hh b/src/hybrid/actor.hh similarity index 100% rename from src/ross_files/actor.hh rename to src/hybrid/actor.hh diff --git a/src/ross_files/generateInfrastructure.py b/src/hybrid/generateInfrastructure.py similarity index 100% rename from src/ross_files/generateInfrastructure.py rename to src/hybrid/generateInfrastructure.py diff --git a/src/ross_files/hybrid.cpp b/src/hybrid/hybrid.cpp similarity index 99% rename from src/ross_files/hybrid.cpp rename to src/hybrid/hybrid.cpp index f67582c..e07c402 100644 --- a/src/ross_files/hybrid.cpp +++ b/src/hybrid/hybrid.cpp @@ -26,8 +26,8 @@ #include "utils/to_string.hpp" #include "reaction_network/network.hpp" #include "hybrid.hpp" -#include "ross_files/simtime.hh" -#include "ross_files/Funnel.hh" +#include "hybrid/simtime.hh" +#include "hybrid/Funnel.hh" #define OPTIONS "hi:s:t:" static const struct option longopts[] = { diff --git a/src/ross_files/hybrid.hpp b/src/hybrid/hybrid.hpp similarity index 100% rename from src/ross_files/hybrid.hpp rename to src/hybrid/hybrid.hpp diff --git a/src/ross_files/simtime.hh b/src/hybrid/simtime.hh similarity index 100% rename from src/ross_files/simtime.hh rename to src/hybrid/simtime.hh From 762f10994a434342f272fca333231be30a6ef412 Mon Sep 17 00:00:00 2001 From: Konstantia Georgouli Date: Mon, 13 Sep 2021 09:59:15 -0700 Subject: [PATCH 08/14] Put in capitals forward, backward, commit and checkpoint in hybrid folder files to avoid any conflicts --- src/hybrid/Checkpoint.hh | 12 +-- src/hybrid/Funnel.hh | 64 +++++++----- src/hybrid/actor.hh | 50 ++++----- src/hybrid/generateInfrastructure.py | 36 +++---- src/hybrid/hybrid.cpp | 149 +++++++++++++++++++++------ 5 files changed, 208 insertions(+), 103 deletions(-) diff --git a/src/hybrid/Checkpoint.hh b/src/hybrid/Checkpoint.hh index 092ad0b..cee1f6f 100644 --- a/src/hybrid/Checkpoint.hh +++ b/src/hybrid/Checkpoint.hh @@ -7,7 +7,7 @@ // This could be done any number of ways faster/better than a list, // but keep it for now. We can optimize this later, so long as we // have the following interface. -class CheckpointQueue { +class CHECKPOINTQueue { public: void push_back(std::string newString) { queue.push_back(newString); } void pop_back() { queue.pop_back(); } @@ -18,12 +18,12 @@ class CheckpointQueue { std::list queue; }; -class InputCheckpoint { +class InputCHECKPOINT { public: std::string str() { return sss.str(); } template - InputCheckpoint& operator&(const TTT& ttt) { + InputCHECKPOINT& operator&(const TTT& ttt) { sss.write(reinterpret_cast(&ttt), sizeof(TTT)); return *this; } @@ -31,13 +31,13 @@ class InputCheckpoint { std::stringstream sss; }; -class OutputCheckpoint { +class OutputCHECKPOINT { public: - OutputCheckpoint(std::string newString) { + OutputCHECKPOINT(std::string newString) { sss.str(newString); } template - OutputCheckpoint& operator&(TTT& ttt) { + OutputCHECKPOINT& operator&(TTT& ttt) { sss.read(reinterpret_cast(&ttt), sizeof(TTT)); return *this; } diff --git a/src/hybrid/Funnel.hh b/src/hybrid/Funnel.hh index 3283f7b..675c165 100644 --- a/src/hybrid/Funnel.hh +++ b/src/hybrid/Funnel.hh @@ -1,14 +1,32 @@ #pragma once #include "actor.hh" +// #include "../reaction_network/network.hpp" +// #include "utils/graph_factory.hpp" +// #include #include #include +#include + +#define MAX_TIME 100 //FIXME! + +// using graph_t = boost::adjacency_list; +/// The type of BGL vertex descriptor for graph_t +// using v_desc_t = boost::graph_traits::vertex_descriptor; +// using v_desc_t = boost::graph_traits::vertex_descriptor; + -#define MAX_TIME 100 //FIXME! typedef std::size_t SpeciesTag; typedef std::size_t ReactionTag; typedef double Real; typedef std::size_t SpeciesValue; + +// typedef reaction_rate_t ( * rate_function_pointer)(const std::vector&); +// typedef std::function&)> RateFunction; + +// typedef std::function&)> f_rate typedef Real (*RateFunction)(const std::vector& species); struct NoopMsg { @@ -36,23 +54,23 @@ class NextReactionMethodCrucible : public LP { ACTOR_DEF(); SLOT(firedMyself, NoopMsg); - void forward(firedMyself)(Time tnow, NoopMsg& msg, tw_bf*) { + void FORWARD(firedMyself)(Time tnow, NoopMsg& msg, tw_bf*) { ReactionMsg r; r.reaction = _rxntag; fire(0, r); _countdown = newCountdown(); fireNextReaction(tnow); } - void commit(firedMyself)(Time tnow, NoopMsg& msg, tw_bf*) { + void COMMIT(firedMyself)(Time tnow, NoopMsg& msg, tw_bf*) { //std::cout << tnow << ": Fired Reaction " << _rxntag << "!\n"; } template - void checkpoint(firedMyself)(Checkpointer& ar, Time tnow, NoopMsg& msg, tw_bf*) { + void CHECKPOINT(firedMyself)(Checkpointer& ar, Time tnow, NoopMsg& msg, tw_bf*) { ar & _countdown & _messageOutstanding & _lastUpdatedTime & _simpleRand; } SLOT(updatedRate, RateMsg); - void forward(updatedRate)(Time tnow, RateMsg& msg, tw_bf*) { + void FORWARD(updatedRate)(Time tnow, RateMsg& msg, tw_bf*) { if (_rate != 0) { EVENT_CANCEL(_messageOutstanding); _countdown -= countdownChange(tnow); @@ -63,7 +81,7 @@ class NextReactionMethodCrucible : public LP { } } template - void checkpoint(updatedRate)(Checkpointer& ar, Time tnow, RateMsg& msg, tw_bf*) { + void CHECKPOINT(updatedRate)(Checkpointer& ar, Time tnow, RateMsg& msg, tw_bf*) { ar & _messageOutstanding & _countdown & _rate & _lastUpdatedTime & _simpleRand; } @@ -145,7 +163,7 @@ class Funnel : public LP { SLOT(firedReaction, ReactionMsg); - void forward(firedReaction)(Time tnow, ReactionMsg& msg, tw_bf* twbf) { + void FORWARD(firedReaction)(Time tnow, ReactionMsg& msg, tw_bf* twbf) { //update species for (const auto& kv : stoic[msg.reaction]) { auto iter = _indexFromSpecTag.find(kv.first); @@ -160,28 +178,28 @@ class Funnel : public LP { setupFutureEvents(tnow, true); } } - void backward(firedReaction)(Time tnow, ReactionMsg& msg, tw_bf*) { + void BACKWARD(firedReaction)(Time tnow, ReactionMsg& msg, tw_bf*) { } - void commit(firedReaction)(Time tnow, ReactionMsg& msg, tw_bf*) { + void COMMIT(firedReaction)(Time tnow, ReactionMsg& msg, tw_bf*) { //std::cout << tnow << ": NRM " << msg.reaction << "->" << _rxntag << std::endl; //printSpecies(tnow); } template - void checkpoint(firedReaction)(Checkpointer& ar, Time tnow, ReactionMsg& msg, tw_bf*) { + void CHECKPOINT(firedReaction)(Checkpointer& ar, Time tnow, ReactionMsg& msg, tw_bf*) { } SLOT(changedRate, FunnelRateExchangeMsg); - void forward(changedRate)(Time tnow, FunnelRateExchangeMsg& msg, tw_bf* twbf) { + void FORWARD(changedRate)(Time tnow, FunnelRateExchangeMsg& msg, tw_bf* twbf) { updateSpeciesBasedOnRates(tnow); updateSavedRates(msg.reaction, msg.nrmMode, msg.rate, msg.normalContrib); interruptNextWakeup(tnow); //std::cout << tnow << ": %%% " << ((tnow-_lastRecomputeTime)>=_maxTimestep) << " " << checkRateBounds(tnow) << std::endl; setupFutureEvents(tnow, (tnow-_lastRecomputeTime)>=_maxTimestep || !checkRateBounds(tnow)); } - void backward(changedRate)(Time tnow, FunnelRateExchangeMsg& msg, tw_bf*) { + void BACKWARD(changedRate)(Time tnow, FunnelRateExchangeMsg& msg, tw_bf*) { } - void commit(changedRate)(Time tnow, FunnelRateExchangeMsg& msg, tw_bf*) { + void COMMIT(changedRate)(Time tnow, FunnelRateExchangeMsg& msg, tw_bf*) { //std::cout << tnow << ": Changed rate for " << msg.reaction << "->" // << _rxntag << ": " // << (msg.nrmMode ? "NRM " : "SDE ") @@ -190,27 +208,27 @@ class Funnel : public LP { //printSpecies(tnow); } template - void checkpoint(changedRate)(Checkpointer& ar, Time tnow, FunnelRateExchangeMsg& msg, tw_bf*) { + void CHECKPOINT(changedRate)(Checkpointer& ar, Time tnow, FunnelRateExchangeMsg& msg, tw_bf*) { } SLOT(wakeup, NoopMsg); - void forward(wakeup)(Time tnow, NoopMsg& msg, tw_bf* twbf) { + void FORWARD(wakeup)(Time tnow, NoopMsg& msg, tw_bf* twbf) { // This gets called when we're sure that we're going to go // out-of-bounds due to our tau condition, based on what we // currently know about the SDE rates. updateSpeciesBasedOnRates(tnow); setupFutureEvents(tnow, (tnow-_lastRecomputeTime)>=_maxTimestep || !checkRateBounds(tnow)); } - void commit(wakeup)(Time tnow, NoopMsg& msg, tw_bf*) { + void COMMIT(wakeup)(Time tnow, NoopMsg& msg, tw_bf*) { //std::cout << tnow << ": Woke Reaction " << _rxntag << "!\n"; //printSpecies(tnow); } template - void checkpoint(wakeup)(Checkpointer& ar, Time tnow, NoopMsg& msg, tw_bf*) { + void CHECKPOINT(wakeup)(Checkpointer& ar, Time tnow, NoopMsg& msg, tw_bf*) { } SLOT(changedMode, SwitchModeMsg); - void forward(changedMode)(Time tnow, SwitchModeMsg& msg, tw_bf*) { + void FORWARD(changedMode)(Time tnow, SwitchModeMsg& msg, tw_bf*) { updateSpeciesBasedOnRates(tnow); int irxn = _indexFromRxnTag[msg.reaction]; updateSavedRates(msg.reaction, msg.nrmMode, _rxnRate[irxn], _normalContrib[irxn]); @@ -220,17 +238,17 @@ class Funnel : public LP { setupFutureEvents(tnow, true); } } - void commit(changedMode)(Time tnow, SwitchModeMsg& msg, tw_bf*) { + void COMMIT(changedMode)(Time tnow, SwitchModeMsg& msg, tw_bf*) { //std::cout << tnow << ": Changed mode " << msg.reaction << "->" << _rxntag << ": " << msg.nrmMode << std::endl; } template - void checkpoint(changedMode)(Checkpointer& ar, Time tnow, SwitchModeMsg& msg, tw_bf*) { + void CHECKPOINT(changedMode)(Checkpointer& ar, Time tnow, SwitchModeMsg& msg, tw_bf*) { } SLOT(printData, NoopMsg); - void forward(printData)(Time tnow, NoopMsg& msg, tw_bf*) { + void FORWARD(printData)(Time tnow, NoopMsg& msg, tw_bf*) { } - void commit(printData)(Time tnow, NoopMsg& msg, tw_bf*) { + void COMMIT(printData)(Time tnow, NoopMsg& msg, tw_bf*) { printSpecies(tnow); } @@ -549,7 +567,7 @@ class Heartbeat : public LP { void setInterval(double interval) { _interval = interval; } SLOT(activated, NoopMsg); - void forward(activated)(Time tnow, NoopMsg& msg, tw_bf*) { + void FORWARD(activated)(Time tnow, NoopMsg& msg, tw_bf*) { activate(_interval, msg); } diff --git a/src/hybrid/actor.hh b/src/hybrid/actor.hh index 5dadc26..082976b 100644 --- a/src/hybrid/actor.hh +++ b/src/hybrid/actor.hh @@ -29,59 +29,59 @@ class BaseLP { tw_lpid getID() const; }; -extern CheckpointQueue checkpointQueue; +extern CHECKPOINTQueue checkpointQueue; template class LP : public BaseLP { public: typedef void (Dependent::*MethodPtr)(Time tnow, void*, tw_bf*); static void init_lp(void *state, tw_lp *me); - static void forward_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp); - static void backward_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp); - static void commit_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp); + static void FORWARD_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp); + static void BACKWARD_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp); + static void COMMIT_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp); static void final_lp(void *state, tw_lp *me); - static MethodPtr forwardDispatch[]; - static MethodPtr backwardDispatch[]; - static MethodPtr commitDispatch[]; + static MethodPtr FORWARDDispatch[]; + static MethodPtr BACKWARDDispatch[]; + static MethodPtr COMMITDispatch[]; template - void forwardWithCheckpoint(Time tnow, typename SlotType::MsgType& msg, tw_bf* twbf) { - InputCheckpoint ic; - static_cast(this)->Dependent::template checkpointFull(ic, tnow, msg, twbf); - static_cast(this)->Dependent::template forwardFull(tnow, msg, twbf); + void FORWARDWithCHECKPOINT(Time tnow, typename SlotType::MsgType& msg, tw_bf* twbf) { + InputCHECKPOINT ic; + static_cast(this)->Dependent::template CHECKPOINTFull(ic, tnow, msg, twbf); + static_cast(this)->Dependent::template FORWARDFull(tnow, msg, twbf); checkpointQueue.push_back(ic.str()); } template - void backwardWithCheckpoint(Time tnow, typename SlotType::MsgType& msg, tw_bf* twbf) { - static_cast(this)->Dependent::template backwardFull(tnow, msg, twbf); - OutputCheckpoint oc(checkpointQueue.back()); - static_cast(this)->Dependent::template checkpointFull(oc, tnow, msg, twbf); + void BACKWARDWithCHECKPOINT(Time tnow, typename SlotType::MsgType& msg, tw_bf* twbf) { + static_cast(this)->Dependent::template BACKWARDFull(tnow, msg, twbf); + OutputCHECKPOINT oc(checkpointQueue.back()); + static_cast(this)->Dependent::template CHECKPOINTFull(oc, tnow, msg, twbf); checkpointQueue.pop_back(); } template - void commitWithCheckpoint(Time tnow, typename SlotType::MsgType& msg, tw_bf* twbf) { - static_cast(this)->Dependent::template commitFull(tnow, msg, twbf); + void COMMITWithCHECKPOINT(Time tnow, typename SlotType::MsgType& msg, tw_bf* twbf) { + static_cast(this)->Dependent::template COMMITFull(tnow, msg, twbf); checkpointQueue.pop_front(); } }; #define ACTOR_DEF() \ template \ - void forwardFull(Time tnow, typename SlotType::MsgType& msg, tw_bf*); \ + void FORWARDFull(Time tnow, typename SlotType::MsgType& msg, tw_bf*); \ template \ - void backwardFull(Time tnow, typename SlotType::MsgType& msg, tw_bf*) {} \ + void BACKWARDFull(Time tnow, typename SlotType::MsgType& msg, tw_bf*) {} \ template \ - void commitFull(Time tnow, typename SlotType::MsgType& msg, tw_bf*) {} \ + void COMMITFull(Time tnow, typename SlotType::MsgType& msg, tw_bf*) {} \ template \ - inline void checkpointFull(Archiver& ar, Time tnow, typename SlotType::MsgType& msg, tw_bf*); + inline void CHECKPOINTFull(Archiver& ar, Time tnow, typename SlotType::MsgType& msg, tw_bf*); #define TAGIFY(SlotType) SlotType::tag -#define forward(SlotType) forward_##SlotType -#define backward(SlotType) backward_##SlotType -#define commit(SlotType) commit_##SlotType -#define checkpoint(SlotType) checkpoint_##SlotType +#define FORWARD(SlotType) FORWARD_##SlotType +#define BACKWARD(SlotType) BACKWARD_##SlotType +#define COMMIT(SlotType) COMMIT_##SlotType +#define CHECKPOINT(SlotType) CHECKPOINT_##SlotType #define SLOT(name, msgType) \ class name { \ diff --git a/src/hybrid/generateInfrastructure.py b/src/hybrid/generateInfrastructure.py index 6a968f4..e5cbc1c 100755 --- a/src/hybrid/generateInfrastructure.py +++ b/src/hybrid/generateInfrastructure.py @@ -42,7 +42,7 @@ def main(): slotName = mmm.group(1) msgName = mmm.group(2) actors[actorName].addSlot(Slot(slotName,msgName)) - for eventType in ['forward','backward','commit','checkpoint']: + for eventType in ['FORWARD','BACKWARD','COMMIT','CHECKPOINT']: mmm = re.search("void "+eventType+r'\(([A-Za-z0-9_]+)\)',line) if mmm: slotName = mmm.group(1) @@ -75,7 +75,7 @@ def main(): return sizeof(EventData); } -CheckpointQueue checkpointQueue; +CHECKPOINTQueue checkpointQueue; template PdesEvent my_event_send(tw_lpid dest, Time trecv, tw_lp* twlp, Tag tag, TTT& msg) @@ -100,23 +100,23 @@ def main(): "slot" : slot.name, "msg" : slot.msg, } - for eventType in "forward","backward","commit": + for eventType in "FORWARD","BACKWARD","COMMIT": ttt['event'] = eventType if slot.hasMode(eventType): print("""template <> inline void %(actor)s::%(event)sFull<%(actor)s::%(slot)s>(Time tnow, %(msg)s& msg, tw_bf* twbf) { %(event)s(%(slot)s)(tnow,msg,twbf); }""" % ttt) - if slot.hasMode('checkpoint'): + if slot.hasMode('CHECKPOINT'): print("""template <> -inline void %(actor)s::checkpointFull<%(actor)s::%(slot)s>(InputCheckpoint& ar, Time tnow, %(msg)s& msg, tw_bf* twbf) { checkpoint(%(slot)s)(ar,tnow,msg,twbf); } +inline void %(actor)s::CHECKPOINTFull<%(actor)s::%(slot)s>(InputCHECKPOINT& ar, Time tnow, %(msg)s& msg, tw_bf* twbf) { CHECKPOINT(%(slot)s)(ar,tnow,msg,twbf); } template <> -inline void %(actor)s::checkpointFull<%(actor)s::%(slot)s>(OutputCheckpoint& ar, Time tnow, %(msg)s& msg, tw_bf* twbf) { checkpoint(%(slot)s)(ar,tnow,msg,twbf); } +inline void %(actor)s::CHECKPOINTFull<%(actor)s::%(slot)s>(OutputCHECKPOINT& ar, Time tnow, %(msg)s& msg, tw_bf* twbf) { CHECKPOINT(%(slot)s)(ar,tnow,msg,twbf); } """ % ttt) for actor in actors.values(): for (ii,slot) in enumerate(actor.getSlots()): print("Tag %s::%s::tag = %d;" % (actor.name,slot.name,ii)) for actor in actors.values(): - for eventType in ["forward","backward","commit"]: + for eventType in ["FORWARD","BACKWARD","COMMIT"]: print("template<> %s::MethodPtr LP<%s>::%sDispatch[] = {" % (actor.name, actor.name,eventType)) for slot in actor.getSlots(): ttt = { @@ -124,8 +124,8 @@ def main(): "slot" : slot.name, "event" : eventType, } - if slot.hasMode('checkpoint'): - method = '%(actor)s::%(event)sWithCheckpoint<%(actor)s::%(slot)s>' % ttt + if slot.hasMode('CHECKPOINT'): + method = '%(actor)s::%(event)sWithCHECKPOINT<%(actor)s::%(slot)s>' % ttt else: method = '%(actor)s::%(event)sFull<%(actor)s::%(slot)s>' % ttt ttt['method'] = method @@ -143,27 +143,27 @@ def main(): } template -void LP::forward_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp) { +void LP::FORWARD_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp) { Time tnow = tw_now(thislp); EventData *e = reinterpret_cast(evt); Dependent* derivedObj = reinterpret_cast(state); - (derivedObj->*forwardDispatch[e->tag])(tnow, &e->msg, bf); + (derivedObj->*FORWARDDispatch[e->tag])(tnow, &e->msg, bf); } template -void LP::backward_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp) { +void LP::BACKWARD_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp) { Time tnow = tw_now(thislp); EventData *e = reinterpret_cast(evt); Dependent* derivedObj = reinterpret_cast(state); - (derivedObj->*backwardDispatch[e->tag])(tnow, &e->msg, bf); + (derivedObj->*BACKWARDDispatch[e->tag])(tnow, &e->msg, bf); } template -void LP::commit_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp) { +void LP::COMMIT_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp) { Time tnow = tw_now(thislp); EventData *e = reinterpret_cast(evt); Dependent* derivedObj = reinterpret_cast(state); - (derivedObj->*commitDispatch[e->tag])(tnow, &e->msg, bf); + (derivedObj->*COMMITDispatch[e->tag])(tnow, &e->msg, bf); } template @@ -175,9 +175,9 @@ def main(): """) for actor in actors.values(): print(""" -template void LP<%(actor)s>::forward_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp); -template void LP<%(actor)s>::backward_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp); -template void LP<%(actor)s>::commit_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp); +template void LP<%(actor)s>::FORWARD_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp); +template void LP<%(actor)s>::BACKWARD_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp); +template void LP<%(actor)s>::COMMIT_event(void *state, tw_bf *bf, void *evt, tw_lp *thislp); template void LP<%(actor)s>::final_lp(void *state, tw_lp *me); template void LP<%(actor)s>::init_lp(void *state, tw_lp *me); diff --git a/src/hybrid/hybrid.cpp b/src/hybrid/hybrid.cpp index e07c402..32b0374 100644 --- a/src/hybrid/hybrid.cpp +++ b/src/hybrid/hybrid.cpp @@ -21,10 +21,14 @@ #include #include #include +// #include "bgl.hpp" #include "utils/write_graphviz.hpp" #include "utils/timer.hpp" #include "utils/to_string.hpp" #include "reaction_network/network.hpp" +#include "reaction_network/species.hpp" +#include "reaction_network/reaction.hpp" +#include "reaction_network/vertex.hpp" #include "hybrid.hpp" #include "hybrid/simtime.hh" #include "hybrid/Funnel.hh" @@ -38,6 +42,30 @@ static const struct option longopts[] = { { 0, 0, 0, 0 }, }; +// using r_prop_t = wcs::Network::r_prop_t; +// /// The type of the BGL graph to represent reaction networks +// using graph_t = boost::adjacency_list< +// wcs::wcs_out_edge_list_t, +// wcs::wcs_vertex_list_t, +// boost::bidirectionalS, +// wcs::Vertex, // vertex property bundle +// wcs::Edge, // edge property bundle +// boost::no_property, +// boost::vecS>; +// /// The type of BGL vertex descriptor for graph_t +// using v_desc_t = boost::graph_traits::vertex_descriptor; +// v_desc_t test; +double interval = 0.2; +int seedBase = 137; +int trust_region = 20000; +int numReactions = 0; +int numSpecies = 0; +const wcs::Network::reaction_list_t* reaction_list; +const wcs::Network::species_list_t* species_list; +std::string filename = ""; +// wcs::Network* rnet; +wcs::Network rnet; + void do_nothing(void *,void *,void *) {} tw_peid ctr_map(tw_lpid gid) { @@ -69,9 +97,9 @@ class RegisterBufferAndOffset { tw_lptype Funnellps[] = { { (init_f) RegisterBufferAndOffset::init_f, (pre_run_f) do_nothing, - (event_f) Funnel::forward_event, - (revent_f) Funnel::backward_event, - (commit_f) Funnel::commit_event, + (event_f) Funnel::FORWARD_event, + (revent_f) Funnel::BACKWARD_event, + (commit_f) Funnel::COMMIT_event, (final_f) Funnel::final_lp, (map_f) ctr_map, //NULL, @@ -87,9 +115,9 @@ int RegisterBufferAndOffset::offset = 0; tw_lptype NRMlps[] = { { (init_f) RegisterBufferAndOffset::init_f, (pre_run_f) do_nothing, - (event_f) NextReactionMethodCrucible::forward_event, - (revent_f) NextReactionMethodCrucible::backward_event, - (commit_f) NextReactionMethodCrucible::commit_event, + (event_f) NextReactionMethodCrucible::FORWARD_event, + (revent_f) NextReactionMethodCrucible::BACKWARD_event, + (commit_f) NextReactionMethodCrucible::COMMIT_event, (final_f) NextReactionMethodCrucible::final_lp, (map_f) ctr_map, //NULL, @@ -107,9 +135,9 @@ int RegisterBufferAndOffset::offset = 0; tw_lptype Heartbeatlps[] = { { (init_f) RegisterBufferAndOffset::init_f, (pre_run_f) do_nothing, - (event_f) Heartbeat::forward_event, - (revent_f) Heartbeat::backward_event, - (commit_f) Heartbeat::commit_event, + (event_f) Heartbeat::FORWARD_event, + (revent_f) Heartbeat::BACKWARD_event, + (commit_f) Heartbeat::COMMIT_event, (final_f) Heartbeat::final_lp, (map_f) ctr_map, //NULL, @@ -141,7 +169,7 @@ enum ReactionEnum { Y_2X_TO_3X, B_X_TO_Y_D, X_TO_E, - numReactions + numReactionss }; enum SpeciesEnum { @@ -151,13 +179,13 @@ enum SpeciesEnum { E_idx, X_idx, Y_idx, - numSpecies + numSpeciess }; // FIXME, comment out, every stoic reference has to check the graph -std::vector> stoic(numReactions); +std::vector> stoic(numReactionss); //FIXME, all of these functions are linked in with dlsym trick @@ -193,8 +221,9 @@ void print_usage(const std::string exec, int code) { std::cerr << "Usage: " << exec << " .xml\n" - " Run the stochastic simulation algorithm (SSA) on a reaction\n" - " network graph (.xml).\n" + " Run the simulation using the hybrid simulation mode which combines \n" + " Stochastic differential equations (SDEs) with stochastic simulation \n" + " algorithm (SSA) on a reaction network graph (.xml).\n" "\n" " OPTIONS:\n" " -h, --help\n" @@ -226,9 +255,9 @@ int main(int argc, char **argv) int c; int rc = EXIT_SUCCESS; std::string outfile; - double interval = 0.2; - int seedBase = 137; - int trust_region = 20000; + // double interval = 0.2; + // int seedBase = 137; + // int trust_region = 20000; while ((c = getopt_long(argc, argv, OPTIONS, longopts, NULL)) != -1) { switch (c) { case 'h': /* --help */ @@ -252,15 +281,22 @@ int main(int argc, char **argv) print_usage (argv[0], 1); } std::string fn(argv[optind]); - std::shared_ptr rnet_ptr = std::make_shared(); - wcs::Network& rnet = *rnet_ptr; + filename = fn; + // std::shared_ptr rnet_ptr = std::make_shared(); + // wcs::Network& rnet = *rnet_ptr; + // wcs::Network rnet; rnet.load(fn); rnet.init(); const wcs::Network::graph_t& g = rnet.graph(); //FIXME, needs to come from graph => DONE - int numReactions = rnet.reaction_list().size(); + + reaction_list = &rnet.reaction_list(); + species_list = &rnet.species_list(); + numReactions = reaction_list->size(); + numSpecies = species_list->size(); + std::cout << "numSpecies: " << numSpecies << std::endl; //printf("Calling tw_init\n"); tw_init(&argc, &argv); @@ -328,11 +364,11 @@ void post_lpinit_setup(tw_pe* pe) { //FIXME, output interval comes from options? - double interval = 0.2; + // double interval = 0.2; heartbeats[0]->setInterval(interval); //FIXME, seed comes from options - int seedBase = 137; + // int seedBase = 137; for (int ireaction=0; ireactionsetTag(ireaction); crucibles[ireaction]->setSeed(seedBase+ireaction); @@ -345,17 +381,68 @@ void post_lpinit_setup(tw_pe* pe) { //FIXME, inititial values are looked up in graph //setup initial conditions SpeciesValue initial_value[numSpecies]; - initial_value[A_idx] = 301*factor; - initial_value[B_idx] = 1806*factor; - initial_value[D_idx] = 0*factor; - initial_value[E_idx] = 0*factor; - initial_value[X_idx] = 1806*factor; - initial_value[Y_idx] = 1806*factor; + // initial_value[A_idx] = 301*factor; + // initial_value[B_idx] = 1806*factor; + // initial_value[D_idx] = 0*factor; + // initial_value[E_idx] = 0*factor; + // initial_value[X_idx] = 1806*factor; + // initial_value[Y_idx] = 1806*factor; + + const wcs::Network::graph_t& g = rnet.graph(); + int i = 0; + for (const auto& sd : rnet.species_list()) { + const auto& sv = g[sd]; // vertex (property) of the species + const auto& sp = sv.property(); // detailed vertex property data + // std::cout << g[sd].get_label() << sp.get_count() << std::endl; + initial_value[i] = sp.get_count(); + i++; + } + for (size_t i=0; i< numSpecies; i++) { + std::cout << initial_value[i] << std::endl; + } + + const wcs::Network::species_list_t& species_list = rnet.species_list(); + const wcs::Network::reaction_list_t& reaction_list = rnet.reaction_list(); + for (const auto& sd : rnet.reaction_list()) { + // const auto& rv = g[sd]; // vertex (property) of the reaction + // const auto& rp = rv.property(); // detailed vertex property data + const auto& rv = g[sd]; // vertex (property) of the reaction + std::cout << rv.get_label() << std::endl; + } //FIXME, numerical parameter, comes from options - int trust_region = 20000; - - + // int trust_region = 20000; + + for (size_t i = 0u; i < reaction_list.size(); ++i) { + std::cout << "reaction" << i << std::endl; + const auto& vd = reaction_list[i]; + // funnels[i]->_rateFunction = vd.ReactionBase::get_calc_rate_fn(); + // reactant species + std::cout << "reactants " << std::endl; + for (const auto ei_in : boost::make_iterator_range(boost::in_edges(vd, g))) { + const auto vd_reactant = boost::source(ei_in, g); + const auto& sv_reactant = g[vd_reactant]; + + const auto& sp_reactant = sv_reactant.property(); + const auto stoichio = g[ei_in].get_stoichiometry_ratio(); + // if (!species_list.get(sv_reactant.get_label())->getConstant()) { + std::cout << sv_reactant.get_label() << stoichio << std::endl; + // funnels[i]->addSpecies(sv_reactant.get_label(), sp_reactant.get_count(), trust_region); + // } + } + std::cout << "products " << std::endl; + // product species + for (const auto ei_out : boost::make_iterator_range(boost::out_edges(vd, g))) { + const auto vd_product = boost::target(ei_out, g); + const auto& sv_product = g[vd_product]; + + const auto& sp_product = sv_product.property(); + const auto stoichio = g[ei_out].get_stoichiometry_ratio(); + std::cout << sv_product.get_label() << stoichio << std::endl; + // funnels[i]->addSpecies(sv_product.get_label(), sp_product.get_count(), trust_region); + } + } + //FIXME, run addSpecies for each dep species, add the rate function. //stoic[A_TO_X][A_idx] = -1; From 6a90d6c5d993da395bd2206d145e52516d0ea18c Mon Sep 17 00:00:00 2001 From: Konstantia Georgouli Date: Wed, 15 Sep 2021 11:11:08 -0700 Subject: [PATCH 09/14] Implement hybrid simulation mode --- src/hybrid/Funnel.hh | 72 ++--- src/hybrid/hybrid.cpp | 351 ++++++++++++++++--------- src/reaction_network/reaction_base.cpp | 8 + src/reaction_network/reaction_base.hpp | 4 + 4 files changed, 273 insertions(+), 162 deletions(-) diff --git a/src/hybrid/Funnel.hh b/src/hybrid/Funnel.hh index 675c165..99fc29c 100644 --- a/src/hybrid/Funnel.hh +++ b/src/hybrid/Funnel.hh @@ -1,33 +1,27 @@ #pragma once #include "actor.hh" -// #include "../reaction_network/network.hpp" +#include "reaction_network/network.hpp" // #include "utils/graph_factory.hpp" // #include #include #include #include +#include "wcs_types.hpp" #define MAX_TIME 100 //FIXME! -// using graph_t = boost::adjacency_list; -/// The type of BGL vertex descriptor for graph_t -// using v_desc_t = boost::graph_traits::vertex_descriptor; -// using v_desc_t = boost::graph_traits::vertex_descriptor; - - -typedef std::size_t SpeciesTag; +// typedef std::size_t SpeciesTag; +typedef wcs::Network::v_desc_t SpeciesTag; typedef std::size_t ReactionTag; -typedef double Real; +// typedef wcs::Network::v_desc_t ReactionTag; +typedef wcs::reaction_rate_t Real; // replace double with reaction_rate_t typedef std::size_t SpeciesValue; -// typedef reaction_rate_t ( * rate_function_pointer)(const std::vector&); -// typedef std::function&)> RateFunction; -// typedef std::function&)> f_rate -typedef Real (*RateFunction)(const std::vector& species); +typedef std::function&)> RateFunction; + +// typedef Real (*RateFunction)(const std::vector& species); struct NoopMsg { }; @@ -200,12 +194,12 @@ class Funnel : public LP { void BACKWARD(changedRate)(Time tnow, FunnelRateExchangeMsg& msg, tw_bf*) { } void COMMIT(changedRate)(Time tnow, FunnelRateExchangeMsg& msg, tw_bf*) { - //std::cout << tnow << ": Changed rate for " << msg.reaction << "->" + // std::cout << tnow << ": Changed rate for " << msg.reaction << "->" // << _rxntag << ": " // << (msg.nrmMode ? "NRM " : "SDE ") // << msg.rate << " " // << msg.normalContrib << std::endl; - //printSpecies(tnow); + printSpecies(tnow); } template void CHECKPOINT(changedRate)(Checkpointer& ar, Time tnow, FunnelRateExchangeMsg& msg, tw_bf*) { @@ -256,7 +250,7 @@ class Funnel : public LP { SIGNAL(updateRateFunnel, FunnelRateExchangeMsg); SIGNAL(updateModeFunnel, SwitchModeMsg); - void addSpecies(int tag, SpeciesValue initialValue, int trustRegion) { + void addSpecies(SpeciesTag tag, SpeciesValue initialValue, int trustRegion) { int oldSize = _indexFromSpecTag.size(); _indexFromSpecTag[tag] = oldSize; _tagFromSpecIndex.push_back(tag); @@ -267,7 +261,7 @@ class Funnel : public LP { _trustRegion.push_back(trustRegion); } - void setTag(int tag) { + void setTag(ReactionTag tag) { _rxntag = tag; } @@ -310,17 +304,33 @@ class Funnel : public LP { } void printSpecies(Time tnow) { - std::cout << tnow << ": " - << "Reaction " << _rxntag << ": (wakeup " << _wakeupTime << ") (nrm=" << _rxnNrmMode[_thisrxnindex] << " rate=" << _rxnRate[_thisrxnindex] << ")" ; - for (std::size_t ii=0; ii<_tagFromSpecIndex.size(); ii++) { - std::cout << " {" << _tagFromSpecIndex[ii] << ": " - << "(" << _lowerBound[ii] << "<" - << _species[ii] - << "<" << _upperBound[ii] << ") " - << _speciesRate[ii] - << "}"; + // std::cout << tnow << ": " + // << "Reaction " << _rxntag << ": (wakeup " << _wakeupTime << ") (nrm=" << _rxnNrmMode[_thisrxnindex] << " rate=" << _rxnRate[_thisrxnindex] << ")" ; + // for (std::size_t ii=0; ii<_tagFromSpecIndex.size(); ii++) { + // std::cout << " {" << _tagFromSpecIndex[ii] << ": " + // << "(" << _lowerBound[ii] << "<" + // << _species[ii] + // << "<" << _upperBound[ii] << ") " + // << _speciesRate[ii] + // << "}"; + // } + // std::cout << std::endl; + int both_species = 0; + for (int ii=0; ii<_tagFromSpecIndex.size(); ii++) { + if (_tagFromSpecIndex[ii] == 2 || _tagFromSpecIndex[ii] == 4){ + both_species = both_species +1; + } } - std::cout << std::endl; + if (both_species == 2 ) { + std::cout << tnow << " "; + for (int ii=0; ii<_tagFromSpecIndex.size(); ii++) { + if (_tagFromSpecIndex[ii] == 2 || _tagFromSpecIndex[ii] == 4){ + std::cout << _tagFromSpecIndex[ii] << ": " << _species[ii] << " "; + } + } + std::cout << std::endl; + } + } public: @@ -335,7 +345,7 @@ class Funnel : public LP { private: std::unordered_map _indexFromSpecTag; // internal index from tag std::vector _tagFromSpecIndex; // tag from internal index - std::vector _species; //species counts + std::vector _species; //species counts std::vector _speciesRate; //first derivative of the species std::vector _lowerBound; // lb of the tau independence region std::vector _upperBound; // ub of the tau independence region @@ -369,7 +379,7 @@ class Funnel : public LP { //compute gradient std::vector gradient(_species.size()); for (unsigned int ii=0; ii<_species.size(); ii++) { - std::vector modSpecies(_species); + std::vector modSpecies(_species); modSpecies[ii]++; Real gradRate = _rateFunction(modSpecies); gradient[ii] = gradRate-newRate; diff --git a/src/hybrid/hybrid.cpp b/src/hybrid/hybrid.cpp index 32b0374..b3ee948 100644 --- a/src/hybrid/hybrid.cpp +++ b/src/hybrid/hybrid.cpp @@ -33,35 +33,29 @@ #include "hybrid/simtime.hh" #include "hybrid/Funnel.hh" -#define OPTIONS "hi:s:t:" +#define OPTIONS "hi:s:r:t:" static const struct option longopts[] = { {"help", no_argument, 0, 'h'}, {"interval", required_argument, 0, 'i'}, {"seedbase", required_argument, 0, 's'}, - {"trust_region", required_argument, 0, 't'}, - { 0, 0, 0, 0 }, + {"trust_region", required_argument, 0, 'r'}, + {"simulation_time", required_argument, 0, 't'}, + { 0, 0, 0, 0}, }; -// using r_prop_t = wcs::Network::r_prop_t; -// /// The type of the BGL graph to represent reaction networks -// using graph_t = boost::adjacency_list< -// wcs::wcs_out_edge_list_t, -// wcs::wcs_vertex_list_t, -// boost::bidirectionalS, -// wcs::Vertex, // vertex property bundle -// wcs::Edge, // edge property bundle -// boost::no_property, -// boost::vecS>; +using r_prop_t = wcs::Network::r_prop_t; // /// The type of BGL vertex descriptor for graph_t -// using v_desc_t = boost::graph_traits::vertex_descriptor; +using v_desc_t = boost::graph_traits::vertex_descriptor; // v_desc_t test; double interval = 0.2; int seedBase = 137; int trust_region = 20000; +int simulation_time = 100; int numReactions = 0; int numSpecies = 0; const wcs::Network::reaction_list_t* reaction_list; const wcs::Network::species_list_t* species_list; +const LIBSBML_CPP_NAMESPACE::Model* model; std::string filename = ""; // wcs::Network* rnet; wcs::Network rnet; @@ -163,58 +157,60 @@ tw_petype MyPEType = { -//FIXME, use graph ids -enum ReactionEnum { - A_TO_X, - Y_2X_TO_3X, - B_X_TO_Y_D, - X_TO_E, - numReactionss -}; +// //FIXME, use graph ids +// enum ReactionEnum { +// A_TO_X, +// Y_2X_TO_3X, +// B_X_TO_Y_D, +// X_TO_E, +// numReactionss +// }; -enum SpeciesEnum { - A_idx, - B_idx, - D_idx, - E_idx, - X_idx, - Y_idx, - numSpeciess -}; +// enum SpeciesEnum { +// A_idx, +// B_idx, +// D_idx, +// E_idx, +// X_idx, +// Y_idx, +// numSpeciess +// }; // FIXME, comment out, every stoic reference has to check the graph -std::vector> stoic(numReactionss); +// std::vector> stoic(numReactionss); +std::vector> stoic; + //FIXME, all of these functions are linked in with dlsym trick -int factor = 200; -Real k_A_TO_X = 1; -Real k_Y_2X_TO_3X = 2.7574E-06/factor/factor; -Real k_B_X_TO_Y_D = 0.001660538/factor; -Real k_X_TO_E = 1; - -Real function_A_TO_X(const std::vector& species) { - // for (auto specie : species){ - // std::cout << specie << std::endl; - // } - return k_A_TO_X*species[0]; +// int factor = 200; +// Real k_A_TO_X = 1; +// Real k_Y_2X_TO_3X = 2.7574E-06/factor/factor; +// Real k_B_X_TO_Y_D = 0.001660538/factor; +// Real k_X_TO_E = 1; + +// Real function_A_TO_X(const std::vector& species) { +// // for (auto specie : species){ +// // std::cout << specie << std::endl; +// // } +// return k_A_TO_X*species[0]; -} +// } -Real function_Y_2X_TO_3X(const std::vector& species) { - return k_Y_2X_TO_3X*species[0]*species[0]*species[1]; -} +// Real function_Y_2X_TO_3X(const std::vector& species) { +// return k_Y_2X_TO_3X*species[0]*species[0]*species[1]; +// } -Real function_B_X_TO_Y_D(const std::vector& species) { - return k_B_X_TO_Y_D*species[0]*species[1]; -} +// Real function_B_X_TO_Y_D(const std::vector& species) { +// return k_B_X_TO_Y_D*species[0]*species[1]; +// } -Real function_X_TO_E(const std::vector& species) { - return k_X_TO_E*species[0]; -} +// Real function_X_TO_E(const std::vector& species) { +// return k_X_TO_E*species[0]; +// } void print_usage(const std::string exec, int code) @@ -269,9 +265,12 @@ int main(int argc, char **argv) case 's': /* --seedbase */ seedBase = static_cast(atoi(optarg)); break; - case 't': /* --trust_region */ + case 'r': /* --trust_region */ trust_region = static_cast(atoi(optarg)); break; + case 't': /* --simulation_time */ + simulation_time = static_cast(atoi(optarg)); + break; default: print_usage(argv[0], 1); break; @@ -289,6 +288,12 @@ int main(int argc, char **argv) rnet.init(); const wcs::Network::graph_t& g = rnet.graph(); + LIBSBML_CPP_NAMESPACE::SBMLReader reader; + const LIBSBML_CPP_NAMESPACE::SBMLDocument* document + = reader.readSBML(fn); + + model = document->getModel(); + //FIXME, needs to come from graph => DONE @@ -296,7 +301,8 @@ int main(int argc, char **argv) species_list = &rnet.species_list(); numReactions = reaction_list->size(); numSpecies = species_list->size(); - std::cout << "numSpecies: " << numSpecies << std::endl; + // std::vector> stoic(numReactions); + // std::cout << "numSpecies: " << numSpecies << std::endl; //printf("Calling tw_init\n"); tw_init(&argc, &argv); @@ -321,7 +327,8 @@ int main(int argc, char **argv) tw_lp_settype(2*numReactions, &Heartbeatlps[0]); tw_pe_settype(&MyPEType); - // std::cout << "sampling end: " << g_st_sampling_end < 0) { g_tw_ts_end = g_st_sampling_end; // default 100000 } else { @@ -380,7 +387,7 @@ void post_lpinit_setup(tw_pe* pe) { //FIXME, inititial values are looked up in graph //setup initial conditions - SpeciesValue initial_value[numSpecies]; + // SpeciesValue initial_value[numSpecies]; // initial_value[A_idx] = 301*factor; // initial_value[B_idx] = 1806*factor; // initial_value[D_idx] = 0*factor; @@ -390,88 +397,107 @@ void post_lpinit_setup(tw_pe* pe) { const wcs::Network::graph_t& g = rnet.graph(); - int i = 0; - for (const auto& sd : rnet.species_list()) { - const auto& sv = g[sd]; // vertex (property) of the species - const auto& sp = sv.property(); // detailed vertex property data - // std::cout << g[sd].get_label() << sp.get_count() << std::endl; - initial_value[i] = sp.get_count(); - i++; - } - for (size_t i=0; i< numSpecies; i++) { - std::cout << initial_value[i] << std::endl; - } + // int i = 0; + // for (const auto& sd : rnet.species_list()) { + // const auto& sv = g[sd]; // vertex (property) of the species + // const auto& sp = sv.property(); // detailed vertex property data + // std::cout << g[sd].get_label() << ": " << sp.get_count() << std::endl; + // initial_value[i] = sp.get_count(); + // i++; + // } + // for (size_t i=0; i< numSpecies; i++) { + // std::cout << initial_value[i] << std::endl; + // } const wcs::Network::species_list_t& species_list = rnet.species_list(); const wcs::Network::reaction_list_t& reaction_list = rnet.reaction_list(); - for (const auto& sd : rnet.reaction_list()) { - // const auto& rv = g[sd]; // vertex (property) of the reaction - // const auto& rp = rv.property(); // detailed vertex property data - const auto& rv = g[sd]; // vertex (property) of the reaction - std::cout << rv.get_label() << std::endl; - } + const LIBSBML_CPP_NAMESPACE::ListOfReactions* model_reactions + = model->getListOfReactions(); + const LIBSBML_CPP_NAMESPACE::ListOfSpecies* model_species + = model->getListOfSpecies(); + + // for (size_t ireaction = 0u; ireaction < reaction_list.size(); ++ireaction) { + // const auto& vd = reaction_list[ireaction]; + // std::cout << "Reaction" << ireaction << ": " << vd << std::endl; + // } + + // for (size_t ireaction = 0u; ireaction < reaction_list.size(); ++ireaction) { + // const auto& vd = reaction_list[ireaction]; + // crucibles[vd]->setTag(vd); + // crucibles[vd]->setSeed(seedBase+vd); + + // funnels[vd]->setTag(vd); + // funnels[vd]->setSeed(seedBase+numReactions+vd); + // } + + // for (const auto& sd : rnet.reaction_list()) { + // const auto& rv = g[sd]; // vertex (property) of the reaction + // const auto& rp = rv.property(); // detailed vertex property data + // // funnels[0]->_rateFunction = rp.ReactionBase::get_calc_rate_fn(); + // std::cout << rv.get_label() << std::endl; + // } //FIXME, numerical parameter, comes from options // int trust_region = 20000; for (size_t i = 0u; i < reaction_list.size(); ++i) { - std::cout << "reaction" << i << std::endl; + // std::cout << "reaction" << i << std::endl; const auto& vd = reaction_list[i]; - // funnels[i]->_rateFunction = vd.ReactionBase::get_calc_rate_fn(); + + const auto &reaction = *(model_reactions->get(g[vd].get_label())); + const LIBSBML_CPP_NAMESPACE::ListOfSpeciesReferences* reaction_reactants + = reaction.getListOfReactants(); + const LIBSBML_CPP_NAMESPACE::ListOfSpeciesReferences* reaction_modifiers + = reaction.getListOfModifiers(); + + // std::cout << "Reaction label " << g[vd].get_label() << std::endl; + const auto& rp = g[vd].property(); + funnels[i]->_rateFunction = rp.get_calc_rate_fn(); // reactant species - std::cout << "reactants " << std::endl; + // std::cout << "reactants " << std::endl; for (const auto ei_in : boost::make_iterator_range(boost::in_edges(vd, g))) { const auto vd_reactant = boost::source(ei_in, g); const auto& sv_reactant = g[vd_reactant]; const auto& sp_reactant = sv_reactant.property(); const auto stoichio = g[ei_in].get_stoichiometry_ratio(); - // if (!species_list.get(sv_reactant.get_label())->getConstant()) { - std::cout << sv_reactant.get_label() << stoichio << std::endl; - // funnels[i]->addSpecies(sv_reactant.get_label(), sp_reactant.get_count(), trust_region); - // } - } - std::cout << "products " << std::endl; - // product species - for (const auto ei_out : boost::make_iterator_range(boost::out_edges(vd, g))) { - const auto vd_product = boost::target(ei_out, g); - const auto& sv_product = g[vd_product]; - - const auto& sp_product = sv_product.property(); - const auto stoichio = g[ei_out].get_stoichiometry_ratio(); - std::cout << sv_product.get_label() << stoichio << std::endl; - // funnels[i]->addSpecies(sv_product.get_label(), sp_product.get_count(), trust_region); + // Include only reactants without modifiers + if (reaction_reactants->get(sv_reactant.get_label()) != NULL || reaction_modifiers->get(sv_reactant.get_label()) != NULL) { + // std::cout << sv_reactant.get_label() << "vd: " << vd_reactant << std::endl; + funnels[i]->addSpecies(vd_reactant, sp_reactant.get_count(), trust_region); + // std::cout << sv_reactant.get_label() << stoichio << " " << sp_reactant.get_count() << std::endl; + } } } //FIXME, run addSpecies for each dep species, add the rate function. //stoic[A_TO_X][A_idx] = -1; - stoic[A_TO_X][X_idx] = 1; - funnels[A_TO_X]->_rateFunction = &function_A_TO_X; - funnels[A_TO_X]->addSpecies(A_idx, initial_value[A_idx], trust_region); + // stoic[A_TO_X][X_idx] = 1; + // funnels[A_TO_X]->_rateFunction = &function_A_TO_X; + // funnels[A_TO_X]->addSpecies(A_idx, initial_value[A_idx], trust_region); //funnels[A_TO_X]->addSpecies(X_idx, initial_value[X_idx]); - stoic[Y_2X_TO_3X][X_idx] = -2+3; - stoic[Y_2X_TO_3X][Y_idx] = -1; - funnels[Y_2X_TO_3X]->_rateFunction = function_Y_2X_TO_3X; - funnels[Y_2X_TO_3X]->addSpecies(X_idx, initial_value[X_idx], trust_region); - funnels[Y_2X_TO_3X]->addSpecies(Y_idx, initial_value[Y_idx], trust_region); + // stoic[Y_2X_TO_3X][X_idx] = -2+3; + // stoic[Y_2X_TO_3X][Y_idx] = -1; + // funnels[Y_2X_TO_3X]->_rateFunction = function_Y_2X_TO_3X; + // funnels[Y_2X_TO_3X]->addSpecies(X_idx, initial_value[X_idx], trust_region); + // funnels[Y_2X_TO_3X]->addSpecies(Y_idx, initial_value[Y_idx], trust_region); //stoic[B_X_TO_Y_D][B_idx] = -1; - stoic[B_X_TO_Y_D][X_idx] = -1; - stoic[B_X_TO_Y_D][Y_idx] = 1; + // stoic[B_X_TO_Y_D][X_idx] = -1; + // stoic[B_X_TO_Y_D][Y_idx] = 1; //stoic[B_X_TO_Y_D][D_idx] = 1; - funnels[B_X_TO_Y_D]->_rateFunction = function_B_X_TO_Y_D; - funnels[B_X_TO_Y_D]->addSpecies(B_idx, initial_value[B_idx], trust_region); + // funnels[B_X_TO_Y_D]->_rateFunction = function_B_X_TO_Y_D; + // funnels[B_X_TO_Y_D]->addSpecies(B_idx, initial_value[B_idx], trust_region); //funnels[B_X_TO_Y_D]->addSpecies(D_idx, initial_value[D_idx], trust_region); - funnels[B_X_TO_Y_D]->addSpecies(X_idx, initial_value[X_idx], trust_region); + // funnels[B_X_TO_Y_D]->addSpecies(X_idx, initial_value[X_idx], trust_region); //funnels[B_X_TO_Y_D]->addSpecies(Y_idx, initial_value[Y_idx], trust_region); - stoic[X_TO_E][X_idx] = -1; + // stoic[X_TO_E][X_idx] = -1; //stoic[X_TO_E][E_idx] = 1; - funnels[X_TO_E]->_rateFunction = function_X_TO_E; + // funnels[X_TO_E]->_rateFunction = function_X_TO_E; //funnels[X_TO_E]->addSpecies(E_idx, initial_value[E_idx], trust_region); - funnels[X_TO_E]->addSpecies(X_idx, initial_value[X_idx], trust_region); + // funnels[X_TO_E]->addSpecies(X_idx, initial_value[X_idx], trust_region); //Connect up the funnels to the crucibles they own. for( int ii=0; ii speciesForThisRxn; - //replace stoic instance, possibly do better than a n^2 loop because we have a graph - for (auto iter : stoic[ireaction]) { - speciesForThisRxn.push_back(iter.first); - } + for (size_t ireaction = 0u; ireaction < reaction_list.size(); ++ireaction) { funnels[ireaction]->setupRecvFromReaction(ireaction); - for (int jreaction=0; jreactiondependsOnSpecies(speciesTag)) { - CONNECT(*(crucibles[ireaction]),NextReactionMethodCrucible::fire, - *(funnels[jreaction]),Funnel::firedReaction); - if (ireaction != jreaction) { - funnels[ireaction]->setupSendToReaction(jreaction, funnels[jreaction]->getID()); - funnels[jreaction]->setupRecvFromReaction(ireaction); - } - break; + std::unordered_map dependentSpecies; + std::unordered_map::iterator speciesit; + + const auto& vd = reaction_list[ireaction]; + + // std::cout << g[vd].get_label() << std::endl; + // reactant species + for (const auto ei_in : boost::make_iterator_range(boost::in_edges(vd, g))) { + const auto vd_reactant = boost::source(ei_in, g); + const auto& sv_reactant = g[vd_reactant]; + + const auto stoichio = g[ei_in].get_stoichiometry_ratio(); + // Include only reactants without modifiers + if (stoichio > 0) { + // std::cout << sv_reactant.get_label() << " " << stoichio * -1 << std::endl; + dependentSpecies.insert(std::make_pair(vd_reactant, stoichio * -1)); + } + } + // product species + for (const auto ei_out : boost::make_iterator_range(boost::out_edges(vd, g))) { + const auto vd_product = boost::target(ei_out, g); + const auto& sv_product = g[vd_product]; + + const auto stoichio = g[ei_out].get_stoichiometry_ratio(); + speciesit = dependentSpecies.find(vd_product); + if (speciesit == dependentSpecies.cend()) { + // std::cout << sv_product.get_label() << " " << stoichio << std::endl; + dependentSpecies.insert(std::make_pair(vd_product, stoichio)); + } else { + speciesit->second += stoichio; + // std::cout << sv_product.get_label() << " " << speciesit->second << std::endl; + } + } + // std::unordered_map temp = stoic [ireaction]; + // for (auto tempi : temp) { + // std::cout << tempi.first<< " = " << tempi.second << std::endl; + // } + + // for (auto tempi : dependentSpecies) { + // std::cout << "--" << tempi.first<< " = " << tempi.second << std::endl; + // } + // stoic [ireaction] = dependentSpecies; + stoic.push_back(dependentSpecies); + for (size_t jreaction = 0u; jreaction < reaction_list.size(); ++jreaction) { + // const auto& vdj = reaction_list[jreaction]; + for (auto speciesTag : dependentSpecies) { + if (funnels[jreaction]->dependsOnSpecies(speciesTag.first)) { + CONNECT(*(crucibles[ireaction]),NextReactionMethodCrucible::fire, + *(funnels[jreaction]),Funnel::firedReaction); + if (ireaction != jreaction) { + funnels[ireaction]->setupSendToReaction(jreaction, funnels[jreaction]->getID()); + funnels[jreaction]->setupRecvFromReaction(ireaction); } + break; } + } } - } + } + + //FIXME, use graph instead of bruteforcing all possible connections + //connect the reactions to each other. Do this by scanning the reaction products against + //the dependent reactions. + //using O(n^2) loop here, fix with proper bipartite graph later. + // for (int ireaction=0; ireaction speciesForThisRxn; + // std::cout << "ireaction" << ireaction << std::endl; + // //replace stoic instance, possibly do better than a n^2 loop because we have a graph + // for (auto iter : stoic[ireaction]) { + // speciesForThisRxn.push_back(iter.first); + // std::cout << "iter.first " << iter.first << std::endl; + // } + // funnels[ireaction]->setupRecvFromReaction(ireaction); + // for (int jreaction=0; jreactiondependsOnSpecies(speciesTag)) { + // CONNECT(*(crucibles[ireaction]),NextReactionMethodCrucible::fire, + // *(funnels[jreaction]),Funnel::firedReaction); + // if (ireaction != jreaction) { + // funnels[ireaction]->setupSendToReaction(jreaction, funnels[jreaction]->getID()); + // funnels[jreaction]->setupRecvFromReaction(ireaction); + // } + // break; + // } + // } + // } + // } //connect the printing to the heartbeat for output CONNECT((*heartbeats[0]),Heartbeat::activate,(*heartbeats[0]),Heartbeat::activated); diff --git a/src/reaction_network/reaction_base.cpp b/src/reaction_network/reaction_base.cpp index 91b1ab6..57d2f0a 100644 --- a/src/reaction_network/reaction_base.cpp +++ b/src/reaction_network/reaction_base.cpp @@ -125,6 +125,14 @@ void ReactionBase::set_calc_rate_fn( m_calc_rate = calc_rate; } +const std::function< + reaction_rate_t (const std::vector&) + >& ReactionBase::get_calc_rate_fn() const +{ + return m_calc_rate; +} + + reaction_rate_t ReactionBase::calc_rate(std::vector&& params) { params.push_back(m_rate_const); diff --git a/src/reaction_network/reaction_base.hpp b/src/reaction_network/reaction_base.hpp index be799dd..37b475b 100644 --- a/src/reaction_network/reaction_base.hpp +++ b/src/reaction_network/reaction_base.hpp @@ -49,6 +49,10 @@ class ReactionBase : public VertexPropertyBase { void set_calc_rate_fn(const std::function< reaction_rate_t (const std::vector&) >& calc_rate); + const std::function< + reaction_rate_t (const std::vector&) + >& get_calc_rate_fn() const; + /// Overwrite the reaction rate to a given value void set_rate(const reaction_rate_t rate); reaction_rate_t get_rate() const; From 0dc4b18fc75306e7fd2a2885192e6c929c7e2430 Mon Sep 17 00:00:00 2001 From: Konstantia Georgouli Date: Wed, 15 Sep 2021 15:18:59 -0700 Subject: [PATCH 10/14] Clean the code from any useless comments --- src/hybrid/Funnel.hh | 5 +- src/hybrid/hybrid.cpp | 190 +----------------------------------------- 2 files changed, 2 insertions(+), 193 deletions(-) diff --git a/src/hybrid/Funnel.hh b/src/hybrid/Funnel.hh index 99fc29c..f456e13 100644 --- a/src/hybrid/Funnel.hh +++ b/src/hybrid/Funnel.hh @@ -11,17 +11,14 @@ #define MAX_TIME 100 //FIXME! -// typedef std::size_t SpeciesTag; typedef wcs::Network::v_desc_t SpeciesTag; typedef std::size_t ReactionTag; -// typedef wcs::Network::v_desc_t ReactionTag; -typedef wcs::reaction_rate_t Real; // replace double with reaction_rate_t +typedef wcs::reaction_rate_t Real; typedef std::size_t SpeciesValue; typedef std::function&)> RateFunction; -// typedef Real (*RateFunction)(const std::vector& species); struct NoopMsg { }; diff --git a/src/hybrid/hybrid.cpp b/src/hybrid/hybrid.cpp index b3ee948..4d31304 100644 --- a/src/hybrid/hybrid.cpp +++ b/src/hybrid/hybrid.cpp @@ -21,7 +21,6 @@ #include #include #include -// #include "bgl.hpp" #include "utils/write_graphviz.hpp" #include "utils/timer.hpp" #include "utils/to_string.hpp" @@ -46,7 +45,6 @@ static const struct option longopts[] = { using r_prop_t = wcs::Network::r_prop_t; // /// The type of BGL vertex descriptor for graph_t using v_desc_t = boost::graph_traits::vertex_descriptor; -// v_desc_t test; double interval = 0.2; int seedBase = 137; int trust_region = 20000; @@ -57,7 +55,6 @@ const wcs::Network::reaction_list_t* reaction_list; const wcs::Network::species_list_t* species_list; const LIBSBML_CPP_NAMESPACE::Model* model; std::string filename = ""; -// wcs::Network* rnet; wcs::Network rnet; void do_nothing(void *,void *,void *) {} @@ -156,63 +153,9 @@ tw_petype MyPEType = { }; - -// //FIXME, use graph ids -// enum ReactionEnum { -// A_TO_X, -// Y_2X_TO_3X, -// B_X_TO_Y_D, -// X_TO_E, -// numReactionss -// }; - -// enum SpeciesEnum { -// A_idx, -// B_idx, -// D_idx, -// E_idx, -// X_idx, -// Y_idx, -// numSpeciess -// }; - - - -// FIXME, comment out, every stoic reference has to check the graph -// std::vector> stoic(numReactionss); std::vector> stoic; - -//FIXME, all of these functions are linked in with dlsym trick - -// int factor = 200; -// Real k_A_TO_X = 1; -// Real k_Y_2X_TO_3X = 2.7574E-06/factor/factor; -// Real k_B_X_TO_Y_D = 0.001660538/factor; -// Real k_X_TO_E = 1; - -// Real function_A_TO_X(const std::vector& species) { -// // for (auto specie : species){ -// // std::cout << specie << std::endl; -// // } -// return k_A_TO_X*species[0]; - -// } - -// Real function_Y_2X_TO_3X(const std::vector& species) { -// return k_Y_2X_TO_3X*species[0]*species[0]*species[1]; -// } - -// Real function_B_X_TO_Y_D(const std::vector& species) { -// return k_B_X_TO_Y_D*species[0]*species[1]; -// } - -// Real function_X_TO_E(const std::vector& species) { -// return k_X_TO_E*species[0]; -// } - - void print_usage(const std::string exec, int code) { std::cerr << @@ -251,9 +194,6 @@ int main(int argc, char **argv) int c; int rc = EXIT_SUCCESS; std::string outfile; - // double interval = 0.2; - // int seedBase = 137; - // int trust_region = 20000; while ((c = getopt_long(argc, argv, OPTIONS, longopts, NULL)) != -1) { switch (c) { case 'h': /* --help */ @@ -281,9 +221,6 @@ int main(int argc, char **argv) } std::string fn(argv[optind]); filename = fn; - // std::shared_ptr rnet_ptr = std::make_shared(); - // wcs::Network& rnet = *rnet_ptr; - // wcs::Network rnet; rnet.load(fn); rnet.init(); const wcs::Network::graph_t& g = rnet.graph(); @@ -301,8 +238,6 @@ int main(int argc, char **argv) species_list = &rnet.species_list(); numReactions = reaction_list->size(); numSpecies = species_list->size(); - // std::vector> stoic(numReactions); - // std::cout << "numSpecies: " << numSpecies << std::endl; //printf("Calling tw_init\n"); tw_init(&argc, &argv); @@ -327,7 +262,6 @@ int main(int argc, char **argv) tw_lp_settype(2*numReactions, &Heartbeatlps[0]); tw_pe_settype(&MyPEType); - // std::cout << "sampling end: " << g_st_sampling_end < 0) { g_tw_ts_end = g_st_sampling_end; // default 100000 @@ -370,12 +304,8 @@ void post_lpinit_setup(tw_pe* pe) { Heartbeat** heartbeats = RegisterBufferAndOffset::getBuffer(); - //FIXME, output interval comes from options? - // double interval = 0.2; heartbeats[0]->setInterval(interval); - //FIXME, seed comes from options - // int seedBase = 137; for (int ireaction=0; ireactionsetTag(ireaction); crucibles[ireaction]->setSeed(seedBase+ireaction); @@ -385,29 +315,7 @@ void post_lpinit_setup(tw_pe* pe) { } - //FIXME, inititial values are looked up in graph - //setup initial conditions - // SpeciesValue initial_value[numSpecies]; - // initial_value[A_idx] = 301*factor; - // initial_value[B_idx] = 1806*factor; - // initial_value[D_idx] = 0*factor; - // initial_value[E_idx] = 0*factor; - // initial_value[X_idx] = 1806*factor; - // initial_value[Y_idx] = 1806*factor; - - const wcs::Network::graph_t& g = rnet.graph(); - // int i = 0; - // for (const auto& sd : rnet.species_list()) { - // const auto& sv = g[sd]; // vertex (property) of the species - // const auto& sp = sv.property(); // detailed vertex property data - // std::cout << g[sd].get_label() << ": " << sp.get_count() << std::endl; - // initial_value[i] = sp.get_count(); - // i++; - // } - // for (size_t i=0; i< numSpecies; i++) { - // std::cout << initial_value[i] << std::endl; - // } const wcs::Network::species_list_t& species_list = rnet.species_list(); const wcs::Network::reaction_list_t& reaction_list = rnet.reaction_list(); @@ -416,28 +324,6 @@ void post_lpinit_setup(tw_pe* pe) { const LIBSBML_CPP_NAMESPACE::ListOfSpecies* model_species = model->getListOfSpecies(); - // for (size_t ireaction = 0u; ireaction < reaction_list.size(); ++ireaction) { - // const auto& vd = reaction_list[ireaction]; - // std::cout << "Reaction" << ireaction << ": " << vd << std::endl; - // } - - // for (size_t ireaction = 0u; ireaction < reaction_list.size(); ++ireaction) { - // const auto& vd = reaction_list[ireaction]; - // crucibles[vd]->setTag(vd); - // crucibles[vd]->setSeed(seedBase+vd); - - // funnels[vd]->setTag(vd); - // funnels[vd]->setSeed(seedBase+numReactions+vd); - // } - - // for (const auto& sd : rnet.reaction_list()) { - // const auto& rv = g[sd]; // vertex (property) of the reaction - // const auto& rp = rv.property(); // detailed vertex property data - // // funnels[0]->_rateFunction = rp.ReactionBase::get_calc_rate_fn(); - // std::cout << rv.get_label() << std::endl; - // } - //FIXME, numerical parameter, comes from options - // int trust_region = 20000; for (size_t i = 0u; i < reaction_list.size(); ++i) { // std::cout << "reaction" << i << std::endl; @@ -460,45 +346,13 @@ void post_lpinit_setup(tw_pe* pe) { const auto& sp_reactant = sv_reactant.property(); const auto stoichio = g[ei_in].get_stoichiometry_ratio(); - // Include only reactants without modifiers + // Include only reactants without the constant reactants/modifiers which have been converted into modifiers if (reaction_reactants->get(sv_reactant.get_label()) != NULL || reaction_modifiers->get(sv_reactant.get_label()) != NULL) { - // std::cout << sv_reactant.get_label() << "vd: " << vd_reactant << std::endl; funnels[i]->addSpecies(vd_reactant, sp_reactant.get_count(), trust_region); - // std::cout << sv_reactant.get_label() << stoichio << " " << sp_reactant.get_count() << std::endl; } } } - //FIXME, run addSpecies for each dep species, add the rate function. - - //stoic[A_TO_X][A_idx] = -1; - // stoic[A_TO_X][X_idx] = 1; - // funnels[A_TO_X]->_rateFunction = &function_A_TO_X; - // funnels[A_TO_X]->addSpecies(A_idx, initial_value[A_idx], trust_region); - //funnels[A_TO_X]->addSpecies(X_idx, initial_value[X_idx]); - - // stoic[Y_2X_TO_3X][X_idx] = -2+3; - // stoic[Y_2X_TO_3X][Y_idx] = -1; - // funnels[Y_2X_TO_3X]->_rateFunction = function_Y_2X_TO_3X; - // funnels[Y_2X_TO_3X]->addSpecies(X_idx, initial_value[X_idx], trust_region); - // funnels[Y_2X_TO_3X]->addSpecies(Y_idx, initial_value[Y_idx], trust_region); - - //stoic[B_X_TO_Y_D][B_idx] = -1; - // stoic[B_X_TO_Y_D][X_idx] = -1; - // stoic[B_X_TO_Y_D][Y_idx] = 1; - //stoic[B_X_TO_Y_D][D_idx] = 1; - // funnels[B_X_TO_Y_D]->_rateFunction = function_B_X_TO_Y_D; - // funnels[B_X_TO_Y_D]->addSpecies(B_idx, initial_value[B_idx], trust_region); - //funnels[B_X_TO_Y_D]->addSpecies(D_idx, initial_value[D_idx], trust_region); - // funnels[B_X_TO_Y_D]->addSpecies(X_idx, initial_value[X_idx], trust_region); - //funnels[B_X_TO_Y_D]->addSpecies(Y_idx, initial_value[Y_idx], trust_region); - - // stoic[X_TO_E][X_idx] = -1; - //stoic[X_TO_E][E_idx] = 1; - // funnels[X_TO_E]->_rateFunction = function_X_TO_E; - //funnels[X_TO_E]->addSpecies(E_idx, initial_value[E_idx], trust_region); - // funnels[X_TO_E]->addSpecies(X_idx, initial_value[X_idx], trust_region); - //Connect up the funnels to the crucibles they own. for( int ii=0; ii 0) { - // std::cout << sv_reactant.get_label() << " " << stoichio * -1 << std::endl; dependentSpecies.insert(std::make_pair(vd_reactant, stoichio * -1)); } } @@ -533,22 +385,11 @@ void post_lpinit_setup(tw_pe* pe) { const auto stoichio = g[ei_out].get_stoichiometry_ratio(); speciesit = dependentSpecies.find(vd_product); if (speciesit == dependentSpecies.cend()) { - // std::cout << sv_product.get_label() << " " << stoichio << std::endl; dependentSpecies.insert(std::make_pair(vd_product, stoichio)); } else { speciesit->second += stoichio; - // std::cout << sv_product.get_label() << " " << speciesit->second << std::endl; } } - // std::unordered_map temp = stoic [ireaction]; - // for (auto tempi : temp) { - // std::cout << tempi.first<< " = " << tempi.second << std::endl; - // } - - // for (auto tempi : dependentSpecies) { - // std::cout << "--" << tempi.first<< " = " << tempi.second << std::endl; - // } - // stoic [ireaction] = dependentSpecies; stoic.push_back(dependentSpecies); for (size_t jreaction = 0u; jreaction < reaction_list.size(); ++jreaction) { // const auto& vdj = reaction_list[jreaction]; @@ -566,35 +407,6 @@ void post_lpinit_setup(tw_pe* pe) { } } - //FIXME, use graph instead of bruteforcing all possible connections - //connect the reactions to each other. Do this by scanning the reaction products against - //the dependent reactions. - //using O(n^2) loop here, fix with proper bipartite graph later. - // for (int ireaction=0; ireaction speciesForThisRxn; - // std::cout << "ireaction" << ireaction << std::endl; - // //replace stoic instance, possibly do better than a n^2 loop because we have a graph - // for (auto iter : stoic[ireaction]) { - // speciesForThisRxn.push_back(iter.first); - // std::cout << "iter.first " << iter.first << std::endl; - // } - // funnels[ireaction]->setupRecvFromReaction(ireaction); - // for (int jreaction=0; jreactiondependsOnSpecies(speciesTag)) { - // CONNECT(*(crucibles[ireaction]),NextReactionMethodCrucible::fire, - // *(funnels[jreaction]),Funnel::firedReaction); - // if (ireaction != jreaction) { - // funnels[ireaction]->setupSendToReaction(jreaction, funnels[jreaction]->getID()); - // funnels[jreaction]->setupRecvFromReaction(ireaction); - // } - // break; - // } - // } - // } - // } //connect the printing to the heartbeat for output CONNECT((*heartbeats[0]),Heartbeat::activate,(*heartbeats[0]),Heartbeat::activated); From cd58a345d73de8fbb89bf955ebd0b45f17e26f6c Mon Sep 17 00:00:00 2001 From: Konstantia Georgouli Date: Wed, 15 Sep 2021 15:32:04 -0700 Subject: [PATCH 11/14] Fix some details (printing of arguments and output of Funnel.hh) --- src/hybrid/Funnel.hh | 50 +++++++++++++++++++++---------------------- src/hybrid/hybrid.cpp | 5 ++++- 2 files changed, 29 insertions(+), 26 deletions(-) diff --git a/src/hybrid/Funnel.hh b/src/hybrid/Funnel.hh index f456e13..f3f0bf5 100644 --- a/src/hybrid/Funnel.hh +++ b/src/hybrid/Funnel.hh @@ -301,32 +301,32 @@ class Funnel : public LP { } void printSpecies(Time tnow) { - // std::cout << tnow << ": " - // << "Reaction " << _rxntag << ": (wakeup " << _wakeupTime << ") (nrm=" << _rxnNrmMode[_thisrxnindex] << " rate=" << _rxnRate[_thisrxnindex] << ")" ; - // for (std::size_t ii=0; ii<_tagFromSpecIndex.size(); ii++) { - // std::cout << " {" << _tagFromSpecIndex[ii] << ": " - // << "(" << _lowerBound[ii] << "<" - // << _species[ii] - // << "<" << _upperBound[ii] << ") " - // << _speciesRate[ii] - // << "}"; - // } - // std::cout << std::endl; - int both_species = 0; - for (int ii=0; ii<_tagFromSpecIndex.size(); ii++) { - if (_tagFromSpecIndex[ii] == 2 || _tagFromSpecIndex[ii] == 4){ - both_species = both_species +1; - } - } - if (both_species == 2 ) { - std::cout << tnow << " "; - for (int ii=0; ii<_tagFromSpecIndex.size(); ii++) { - if (_tagFromSpecIndex[ii] == 2 || _tagFromSpecIndex[ii] == 4){ - std::cout << _tagFromSpecIndex[ii] << ": " << _species[ii] << " "; - } - } - std::cout << std::endl; + std::cout << tnow << ": " + << "Reaction " << _rxntag << ": (wakeup " << _wakeupTime << ") (nrm=" << _rxnNrmMode[_thisrxnindex] << " rate=" << _rxnRate[_thisrxnindex] << ")" ; + for (std::size_t ii=0; ii<_tagFromSpecIndex.size(); ii++) { + std::cout << " {" << _tagFromSpecIndex[ii] << ": " + << "(" << _lowerBound[ii] << "<" + << _species[ii] + << "<" << _upperBound[ii] << ") " + << _speciesRate[ii] + << "}"; } + std::cout << std::endl; + // int both_species = 0; + // for (int ii=0; ii<_tagFromSpecIndex.size(); ii++) { + // if (_tagFromSpecIndex[ii] == 2 || _tagFromSpecIndex[ii] == 4){ + // both_species = both_species +1; + // } + // } + // if (both_species == 2 ) { + // std::cout << tnow << " "; + // for (int ii=0; ii<_tagFromSpecIndex.size(); ii++) { + // if (_tagFromSpecIndex[ii] == 2 || _tagFromSpecIndex[ii] == 4){ + // std::cout << _tagFromSpecIndex[ii] << ": " << _species[ii] << " "; + // } + // } + // std::cout << std::endl; + // } } diff --git a/src/hybrid/hybrid.cpp b/src/hybrid/hybrid.cpp index 4d31304..7389669 100644 --- a/src/hybrid/hybrid.cpp +++ b/src/hybrid/hybrid.cpp @@ -174,8 +174,11 @@ void print_usage(const std::string exec, int code) " -s, --seedbase\n" " Specify the seed for random number generator\n" "\n" - " -t, --trust_region\n" + " -r, --trust_region\n" " Specify the trust region\n" + "\n" + " -t, --simulation_time\n" + " Specify the simulation time\n" "\n"; exit(code); } From 1bd0d88a1cd423c3611843892e27168fb20e3ab4 Mon Sep 17 00:00:00 2001 From: Konstantia Georgouli Date: Wed, 22 Sep 2021 10:14:43 -0700 Subject: [PATCH 12/14] Include substance units in the calculations when initial concentration is used --- src/reaction_network/vertex.hpp | 31 ++++++++---- src/utils/generate_cxx_code.cpp | 83 ++++++++++++++++++++++++++------- 2 files changed, 89 insertions(+), 25 deletions(-) diff --git a/src/reaction_network/vertex.hpp b/src/reaction_network/vertex.hpp index 4549cd6..fda8de4 100644 --- a/src/reaction_network/vertex.hpp +++ b/src/reaction_network/vertex.hpp @@ -181,6 +181,7 @@ const LIBSBML_CPP_NAMESPACE::Species& species, const G& g) const LIBSBML_CPP_NAMESPACE::ListOfUnitDefinitions* unit_definition_list = model.getListOfUnitDefinitions(); double compartment_size = compartment_list->get(species.getCompartment())->getSize(); + // Find compartment units std::string compartment_unit ; if (compartment_list->get(species.getCompartment())->isSetUnits()){ compartment_unit = compartment_list->get(species.getCompartment())->getUnits(); @@ -191,19 +192,33 @@ const LIBSBML_CPP_NAMESPACE::Species& species, const G& g) } else if (compartment_list->get(species.getCompartment())->getSpatialDimensions() == 1) { compartment_unit = model.getLengthUnits(); } - const LIBSBML_CPP_NAMESPACE::UnitDefinition* unit_definition = unit_definition_list->get(compartment_unit); - const LIBSBML_CPP_NAMESPACE::ListOfUnits* unit_list - = unit_definition->getListOfUnits(); - unsigned int unitsSize = unit_list->size(); + const LIBSBML_CPP_NAMESPACE::UnitDefinition* unit_definition_comp = unit_definition_list->get(compartment_unit); + const LIBSBML_CPP_NAMESPACE::ListOfUnits* unit_list_comp + = unit_definition_comp->getListOfUnits(); + unsigned int unitsSize_comp = unit_list_comp->size(); double comp_unit = 1.0; - for (unsigned int iu = 0u; iu < unitsSize; iu++) { - const LIBSBML_CPP_NAMESPACE::Unit* unit = unit_list->get(iu); + for (unsigned int iu = 0u; iu < unitsSize_comp; iu++) { + const LIBSBML_CPP_NAMESPACE::Unit* unit = unit_list_comp->get(iu); comp_unit = comp_unit * pow(unit->getMultiplier()*pow(10,unit->getScale()),unit->getExponent()); - } + } + // Find substance units + std::string substance_unit ; + substance_unit = model.getSubstanceUnits(); + + const LIBSBML_CPP_NAMESPACE::UnitDefinition* unit_definition_substance = unit_definition_list->get(substance_unit); + const LIBSBML_CPP_NAMESPACE::ListOfUnits* unit_list_substance + = unit_definition_substance->getListOfUnits(); + unsigned int unitsSize_substance = unit_list_substance->size(); + double sub_unit = 1.0; + for (unsigned int iu = 0u; iu < unitsSize_substance; iu++) { + const LIBSBML_CPP_NAMESPACE::Unit* unit = unit_list_substance->get(iu); + sub_unit = sub_unit * pow(unit->getMultiplier()*pow(10,unit->getScale()),unit->getExponent()); + } + double avog_num = 6.02214e+23; dynamic_cast(m_p.get())-> set_count(static_cast(species.getInitialConcentration() * (6.02214E23 * - compartment_size) *comp_unit)); + compartment_size) * comp_unit * sub_unit)); } } diff --git a/src/utils/generate_cxx_code.cpp b/src/utils/generate_cxx_code.cpp index e1c078c..f29fdea 100644 --- a/src/utils/generate_cxx_code.cpp +++ b/src/utils/generate_cxx_code.cpp @@ -175,7 +175,8 @@ static double count2Concentration_converter ( const LIBSBML_CPP_NAMESPACE::Model& model, const LIBSBML_CPP_NAMESPACE::Species& species, std::string& compartment, - double& comp_unit) + double& comp_unit, + double& sub_unit) { const LIBSBML_CPP_NAMESPACE::ListOfCompartments* compartment_list = model.getListOfCompartments(); @@ -196,17 +197,32 @@ static double count2Concentration_converter ( } else if (compartment_list->get(species.getCompartment())->getSpatialDimensions() == 1) { compartment_unit = model.getLengthUnits(); } - const LIBSBML_CPP_NAMESPACE::UnitDefinition* unit_definition = unit_definition_list->get(compartment_unit); - const LIBSBML_CPP_NAMESPACE::ListOfUnits* unit_list - = unit_definition->getListOfUnits(); - unsigned int unitsSize = unit_list->size(); + const LIBSBML_CPP_NAMESPACE::UnitDefinition* unit_definition_comp = unit_definition_list->get(compartment_unit); + const LIBSBML_CPP_NAMESPACE::ListOfUnits* unit_list_comp + = unit_definition_comp->getListOfUnits(); + unsigned int unitsSize_comp = unit_list_comp->size(); // double comp_unit = 1.0; comp_unit = 1.0; - for (unsigned int iu = 0u; iu < unitsSize; iu++) { - const LIBSBML_CPP_NAMESPACE::Unit* unit = unit_list->get(iu); + for (unsigned int iu = 0u; iu < unitsSize_comp; iu++) { + const LIBSBML_CPP_NAMESPACE::Unit* unit = unit_list_comp->get(iu); comp_unit = comp_unit * pow(unit->getMultiplier()*pow(10,unit->getScale()),unit->getExponent()); - } - return avog_num * compartment_size * comp_unit; + } + + // Find substance units + std::string substance_unit ; + substance_unit = model.getSubstanceUnits(); + + const LIBSBML_CPP_NAMESPACE::UnitDefinition* unit_definition_substance = unit_definition_list->get(substance_unit); + const LIBSBML_CPP_NAMESPACE::ListOfUnits* unit_list_substance + = unit_definition_substance->getListOfUnits(); + unsigned int unitsSize_substance = unit_list_substance->size(); + sub_unit = 1.0; + for (unsigned int iu = 0u; iu < unitsSize_substance; iu++) { + const LIBSBML_CPP_NAMESPACE::Unit* unit = unit_list_substance->get(iu); + sub_unit = sub_unit * pow(unit->getMultiplier()*pow(10,unit->getScale()),unit->getExponent()); + } + + return avog_num * compartment_size * comp_unit * sub_unit; } @@ -1138,6 +1154,9 @@ void generate_cxx_code::print_reaction_rates( typename compartment_t::const_iterator rci; compartment_t reaction_comp; unsigned int num_localparameters = 0u; + double model_comp_unit = 1; + double model_sub_unit = 1; + int model_comp_unit_defined = 0; if (model.getLevel() > 2) { const LIBSBML_CPP_NAMESPACE::ListOfLocalParameters* local_parameter_list = reaction.getKineticLaw()->getListOfLocalParameters(); @@ -1241,8 +1260,18 @@ void generate_cxx_code::print_reaction_rates( if (species_list->get(*itf) != NULL) { if (species_list->get(*itf)->isSetInitialConcentration()) { double comp_unit; + double sub_unit; std::string compart_name; - double convert_factor = count2Concentration_converter(model, *species_list->get(*itf), compart_name, comp_unit); + double convert_factor = count2Concentration_converter(model, *species_list->get(*itf), compart_name, comp_unit, sub_unit); + if (model_comp_unit_defined == 0) { + model_comp_unit = model_comp_unit * comp_unit; + model_sub_unit = model_sub_unit * sub_unit; + model_comp_unit_defined = 1; + } else { + if (model_comp_unit != comp_unit) { + WCS_THROW("Compartments with different units"); + } + } // Add values for the compartment in map rci = reaction_comp.find(compart_name) ; if (rci == reaction_comp.end()) { @@ -1259,8 +1288,18 @@ void generate_cxx_code::print_reaction_rates( if (species_list->get(*itf) != NULL) { if (species_list->get(*itf)->isSetInitialConcentration()) { double comp_unit; + double sub_unit; std::string compart_name; - double convert_factor = count2Concentration_converter(model, *species_list->get(*itf), compart_name, comp_unit); + double convert_factor = count2Concentration_converter(model, *species_list->get(*itf), compart_name, comp_unit, sub_unit); + if (model_comp_unit_defined == 0) { + model_comp_unit = model_comp_unit * comp_unit; + model_sub_unit = model_sub_unit * sub_unit; + model_comp_unit_defined = 1; + } else { + if (model_comp_unit != comp_unit) { + WCS_THROW("Compartments with different units"); + } + } // Add values for the compartment in map rci = reaction_comp.find(compart_name) ; if (rci == reaction_comp.end()) { @@ -1287,8 +1326,18 @@ void generate_cxx_code::print_reaction_rates( if (species_list->get(*it) != NULL) { if (species_list->get(*it)->isSetInitialConcentration()) { double comp_unit; + double sub_unit; std::string compart_name; - double convert_factor = count2Concentration_converter(model, *species_list->get(*it), compart_name, comp_unit); + double convert_factor = count2Concentration_converter(model, *species_list->get(*it), compart_name, comp_unit, sub_unit); + if (model_comp_unit_defined == 0) { + model_comp_unit = model_comp_unit * comp_unit; + model_sub_unit = model_sub_unit * sub_unit; + model_comp_unit_defined = 1; + } else { + if (model_comp_unit != comp_unit) { + WCS_THROW("Compartments with different units"); + } + } // Add values for the compartment in map rci = reaction_comp.find(compart_name) ; if (rci == reaction_comp.end()) { @@ -1437,11 +1486,11 @@ void generate_cxx_code::print_reaction_rates( genfile << " return " << reaction.getIdAttribute() << ";\n"; } else { // Calculate the compartment unit for all the compartments - double comp_unit = 1.0; - for ( auto rci = reaction_comp.begin(); rci != reaction_comp.end(); ++rci ){ - comp_unit = comp_unit * rci->second; - } - genfile << " return " << reaction.getIdAttribute() << " * " << avog_num << " * " << comp_unit << ";\n"; + // double comp_unit = 1.0; + // for ( auto rci = reaction_comp.begin(); rci != reaction_comp.end(); ++rci ){ + // comp_unit = comp_unit * rci->second; + // } + genfile << " return " << reaction.getIdAttribute() << " * " << avog_num << " * " << model_comp_unit << " * " << model_sub_unit << ";\n"; } genfile << "}\n\n"; } From b46f661badee0d2b1ad3d3a8237d6134ce4644db Mon Sep 17 00:00:00 2001 From: Robert Blake Date: Fri, 22 Oct 2021 13:23:09 -0700 Subject: [PATCH 13/14] Reduce 0 time messages. --- src/hybrid/Funnel.hh | 26 +++++++++++++++++++++++++- src/hybrid/hybrid.cpp | 32 ++++++++++++++++++++++++++------ 2 files changed, 51 insertions(+), 7 deletions(-) diff --git a/src/hybrid/Funnel.hh b/src/hybrid/Funnel.hh index f3f0bf5..ebc4684 100644 --- a/src/hybrid/Funnel.hh +++ b/src/hybrid/Funnel.hh @@ -296,10 +296,34 @@ class Funnel : public LP { _lastUpdatedTime = tnow; //recalculate rate - setupFutureEvents(tnow, true); + setupFutureEvents(tnow, false); printSpecies(tnow); } + FunnelRateExchangeMsg initialRate(Time tnow, bool recompute) { + + FunnelRateExchangeMsg msg; + msg.reaction = _rxntag; + + if (recompute) { + msg.rate = recalculateRate(tnow); + } else { + msg.rate = _rxnRate[_thisrxnindex]; + } + Real delta_time; + getNextMode(msg.rate, delta_time, msg.nrmMode); + if (recompute) { + std::normal_distribution<> dis(0,1); + Real unit_normal_random = dis(_simpleRand); + msg.normalContrib = unit_normal_random * std::sqrt(msg.rate); + } + return msg; + } + + void setInitialRate(FunnelRateExchangeMsg& msg) { + updateSavedRates(msg.reaction, msg.nrmMode, msg.rate, msg.normalContrib); + } + void printSpecies(Time tnow) { std::cout << tnow << ": " << "Reaction " << _rxntag << ": (wakeup " << _wakeupTime << ") (nrm=" << _rxnNrmMode[_thisrxnindex] << " rate=" << _rxnRate[_thisrxnindex] << ")" ; diff --git a/src/hybrid/hybrid.cpp b/src/hybrid/hybrid.cpp index 7389669..ae2ced1 100644 --- a/src/hybrid/hybrid.cpp +++ b/src/hybrid/hybrid.cpp @@ -327,6 +327,7 @@ void post_lpinit_setup(tw_pe* pe) { const LIBSBML_CPP_NAMESPACE::ListOfSpecies* model_species = model->getListOfSpecies(); + for (size_t i = 0u; i < reaction_list.size(); ++i) { // std::cout << "reaction" << i << std::endl; @@ -362,6 +363,8 @@ void post_lpinit_setup(tw_pe* pe) { *(crucibles[ii]),NextReactionMethodCrucible::updatedRate); } + std::vector> reactionConnectivity; + for (size_t ireaction = 0u; ireaction < reaction_list.size(); ++ireaction) { funnels[ireaction]->setupRecvFromReaction(ireaction); std::unordered_map dependentSpecies; @@ -394,22 +397,24 @@ void post_lpinit_setup(tw_pe* pe) { } } stoic.push_back(dependentSpecies); + std::set rxnConnectSet; for (size_t jreaction = 0u; jreaction < reaction_list.size(); ++jreaction) { // const auto& vdj = reaction_list[jreaction]; for (auto speciesTag : dependentSpecies) { if (funnels[jreaction]->dependsOnSpecies(speciesTag.first)) { CONNECT(*(crucibles[ireaction]),NextReactionMethodCrucible::fire, *(funnels[jreaction]),Funnel::firedReaction); + rxnConnectSet.insert(jreaction); if (ireaction != jreaction) { funnels[ireaction]->setupSendToReaction(jreaction, funnels[jreaction]->getID()); funnels[jreaction]->setupRecvFromReaction(ireaction); } break; } - } + } } - } - + reactionConnectivity.emplace_back(std::move(rxnConnectSet)); + } //connect the printing to the heartbeat for output CONNECT((*heartbeats[0]),Heartbeat::activate,(*heartbeats[0]),Heartbeat::activated); @@ -418,15 +423,30 @@ void post_lpinit_setup(tw_pe* pe) { *(funnels[ireaction]),Funnel::printData); } + + Time tnow = 0; + //Setup the initial rates: + for (size_t ireaction = 0u; ireaction < reaction_list.size(); ++ireaction) { + auto rate = funnels[ireaction]->initialRate(tnow, true); + for (auto& jreaction : reactionConnectivity[ireaction]) { + funnels[jreaction]->setInitialRate(rate); + } + } + + //Setup the initial modes: + for (size_t ireaction = 0u; ireaction < reaction_list.size(); ++ireaction) { + auto rate = funnels[ireaction]->initialRate(tnow, false); + for (auto& jreaction : reactionConnectivity[ireaction]) { + funnels[jreaction]->setInitialRate(rate); + } + } //setup the initial messages - Time tnow = 0; for (int ireaction=0; ireactionbeginReactions(tnow); - } NoopMsg msg; - heartbeats[0]->activate(tnow, msg); + //heartbeats[0]->activate(tnow, msg); } From 7c8b9e170d4a3b4877bb8753d9fd4d4d115778c4 Mon Sep 17 00:00:00 2001 From: Robert Blake Date: Tue, 14 Dec 2021 10:14:13 -0800 Subject: [PATCH 14/14] Changing types to be more accurate. --- src/hybrid/Funnel.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/hybrid/Funnel.hh b/src/hybrid/Funnel.hh index ebc4684..bde220a 100644 --- a/src/hybrid/Funnel.hh +++ b/src/hybrid/Funnel.hh @@ -247,7 +247,7 @@ class Funnel : public LP { SIGNAL(updateRateFunnel, FunnelRateExchangeMsg); SIGNAL(updateModeFunnel, SwitchModeMsg); - void addSpecies(SpeciesTag tag, SpeciesValue initialValue, int trustRegion) { + void addSpecies(SpeciesTag tag, SpeciesValue initialValue, SpeciesValue trustRegion) { int oldSize = _indexFromSpecTag.size(); _indexFromSpecTag[tag] = oldSize; _tagFromSpecIndex.push_back(tag);