diff --git a/src/reaction_network/network.cpp b/src/reaction_network/network.cpp index 85a91a4..b6d504a 100644 --- a/src/reaction_network/network.cpp +++ b/src/reaction_network/network.cpp @@ -154,9 +154,9 @@ 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; code_generator.generate_code(*model, m_dep_params_f, m_dep_params_nf, m_rate_rules_dep_map); - const std::string library_file = code_generator.compile_code(); using std::operator<<; @@ -192,6 +192,7 @@ void Network::init() m_species.reserve(num_vertices); v_iter_t vi, vi_end; + 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()); @@ -238,6 +239,7 @@ void Network::init() 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))); + } } 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 d394bc4..4549cd6 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 @@ -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. @@ -153,7 +156,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 +166,44 @@ 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)); } } @@ -209,12 +242,31 @@ 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 - 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 +277,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 +330,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(); @@ -317,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 66f935e..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,6 +171,45 @@ generate_cxx_code::get_all_dependencies( return dependencies_no_dupl; } +static double count2Concentration_converter ( + const LIBSBML_CPP_NAMESPACE::Model& model, + const LIBSBML_CPP_NAMESPACE::Species& species, + std::string& compartment, + double& comp_unit) +{ + 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; +} + + 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 = FunctionInterfaceForJIT::SBML_formulaToString_wcs(denominator); + denominators.push_back(new_denominator); } } } @@ -365,15 +406,36 @@ 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() <<": "; + // 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()); + // 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 +1130,55 @@ 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(); + // 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 + = 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()); + // 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 \", \"" @@ -1116,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); @@ -1139,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); } } @@ -1157,6 +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"; + 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); @@ -1199,7 +1343,6 @@ void generate_cxx_code::print_reaction_rates( par_names_nf.push_back(x); } } - // put input parameters into maps dep_params_f.insert(std::make_pair(reaction.getIdAttribute(), par_names)); @@ -1232,8 +1375,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 = FunctionInterfaceForJIT::SBML_formulaToString_wcs(&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 @@ -1289,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"; } } @@ -1674,6 +1827,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