/****************************************************************************** * * Copyright (C) 2020 by * The Salk Institute for Biological Studies * * Use of this source code is governed by an MIT-style * license that can be found in the LICENSE file or at * https://opensource.org/licenses/MIT. * ******************************************************************************/ #include #include "subsystem.h" #include "bng/bng.h" #include "api/component_type.h" #include "api/elementary_molecule_type.h" #include "api/elementary_molecule.h" #include "api/complex.h" #include "api/api_utils.h" using namespace std; namespace MCell { namespace API { void Subsystem::dump() const { std::cout << to_str() << "\n"; } std::shared_ptr Subsystem::get_species_with_id(const species_id_t id) { // not very efficient, we may need some caching/map later for (auto& s: species) { if (s->species_id == id) { return s; } } return std::shared_ptr(nullptr); } void Subsystem::unify_and_register_elementary_molecule_types() { // then go through Species for (std::shared_ptr s: species) { for (size_t i = 0; i < s->elementary_molecules.size(); i++) { // note the reference, we might be changing this spared ptr std::shared_ptr& emt = s->elementary_molecules[i]->elementary_molecule_type; // do we have this exact object already? auto it_ptr_eq = std::find_if( elementary_molecule_types.begin(), elementary_molecule_types.end(), [&emt] (const std::shared_ptr emt_existing) { return emt_existing.get() == emt.get(); // comparing pointers } ); if (it_ptr_eq != elementary_molecule_types.end()) { // same object exists continue; } // do we have an object with the same contents already? auto it_contents_eq = std::find_if( elementary_molecule_types.begin(), elementary_molecule_types.end(), [&emt] (const std::shared_ptr emt_existing) { return *emt_existing == *emt; // comparing contents } ); if (it_contents_eq != elementary_molecule_types.end()) { // same object exists, we must use this one and let shared_ptr ref counting // destroy the previous one (creation of the EMT object to be destroyed // occurs mainly when defining simple species, so the link from species to this emt // is internal and not visible to the user) emt = *it_contents_eq; continue; } // one more option is that some of the elementary molecule types that we encounter are not // fully initialized, for now expecting that the the fully initialized one is already present auto it_name_eq = std::find_if( elementary_molecule_types.begin(), elementary_molecule_types.end(), [&emt] (const std::shared_ptr emt_existing) { return emt_existing->name == emt->name; // comparing contents } ); if (it_name_eq != elementary_molecule_types.end()) { // use the one that is initialized if (!(*it_name_eq)->all_numerical_attributes_are_unset() && emt->all_numerical_attributes_are_unset() ) { // update the pointer used in the Species object emt = *it_name_eq; continue; } if ((*it_name_eq)->all_numerical_attributes_are_unset() && !emt->all_numerical_attributes_are_unset() ) { // update the one in elementary_molecule_types using the Species object *it_name_eq = emt; continue; } // other cases are error and will be reported in add_elementary_molecule_type } // no such object is in the emt list, add it (reports error if such object already exists) add_elementary_molecule_type(emt); } } } void Subsystem::load_bngl_molecule_types_and_reaction_rules( const std::string& file_name, const std::map& parameter_overrides) { BNG::BNGData bng_data; int num_errors = BNG::parse_bngl_file(file_name, bng_data, parameter_overrides); if (num_errors != 0) { throw RuntimeError("Could not parse BNGL file " + file_name + "."); } // now convert everything we parsed into the API classes so that the user can // inspect or manipulate it if needed convert_bng_data_to_subsystem_data(bng_data); } void Subsystem::convert_bng_data_to_subsystem_data(const BNG::BNGData& bng_data) { // elementary molecules for (const BNG::ElemMolType& mt: bng_data.get_elem_mol_types()) { auto res_mt = ElementaryMoleculeType::construct_from_bng_elem_mol_type(bng_data, mt); append_to_vec_canonical_name(elementary_molecule_types, res_mt); } // reaction rules for (const BNG::RxnRule& rr: bng_data.get_rxn_rules()) { convert_reaction_rule(bng_data, rr); } } void Subsystem::convert_reaction_rule(const BNG::BNGData& bng_data, const BNG::RxnRule& bng_rr) { // always a standard rxn if (bng_rr.type != BNG::RxnType::Standard) { throw RuntimeError("Unexpected type of reaction from BNGL file for reaction " + bng_rr.name + "."); } auto res_rr = make_shared(DefaultCtorArgType()); if (bng_rr.name != "") { res_rr->name = bng_rr.name; } // RxnRule res_rr->fwd_rate = bng_rr.base_rate_constant; for (const BNG::Cplx& inst: bng_rr.reactants) { // MCell3R accepts reactants with any orientation res_rr->reactants.push_back( Complex::construct_from_bng_cplx_w_orientation(bng_data, inst, Orientation::ANY)); } for (const BNG::Cplx& inst: bng_rr.products) { // MCell3R always creates products with the orientation up res_rr->products.push_back( Complex::construct_from_bng_cplx_w_orientation(bng_data, inst, Orientation::UP)); } // allow reactions with identical names append_to_vec_canonical_name(reaction_rules, res_rr); } } // namespace API } // namespace MCell