/****************************************************************************** * * 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 #include "mcell4_converter.h" #include "model.h" #include "world.h" #include "release_event.h" #include "clamp_release_event.h" #include "viz_output_event.h" #include "mol_or_rxn_count_event.h" #include "geometry.h" #include "rng.h" #include "isaac64.h" #include "mcell_structs_shared.h" #include "custom_function_call_event.h" #include "bng/bng.h" #include "api/mcell.h" #include "api/bindings.h" #include "api/compartment_utils.h" #include "api/rng_state.h" using namespace std; namespace MCell { namespace API { static viz_mode_t convert_viz_mode(const VizMode m) { switch (m) { case VizMode::ASCII: return ASCII_MODE; case VizMode::CELLBLENDER_V1: return CELLBLENDER_MODE_V1; case VizMode::CELLBLENDER: return CELLBLENDER_MODE_V2; default: throw ValueError("Invalid VizMode value " + to_string((int)m) + "."); } } MCell4Converter::MCell4Converter(Model* model_, World* world_) : model(model_), world(world_), bng_converter(world_->bng_engine.get_data(), world_->bng_engine.get_config()) { } void MCell4Converter::convert_before_init() { assert(model != nullptr); assert(world != nullptr); convert_simulation_setup(); convert_compartments(); convert_elementary_molecule_types(); convert_species(); convert_surface_classes(); convert_rxns(); // find out whether we have a vol vol rxn for all current species (for mcell3 compatibility) // in BNG mode this finds a reaction as well // this is an optimization to tell that we don't need to check the // surroundings of subpartitions in case when there are no vol-vol rxns // WARNING: must be done before geometry objects are converted because wall_subparts_collision_test // depends on it if (!world->get_all_rxns().has_bimol_vol_rxns()) { // the default is true (or read from user) world->config.use_expanded_list = false; } // at this point, we need to create the first (and for now the only) partition // create initial partition with center at 0,0,0 partition_id_t index = world->add_partition(world->config.partition0_llf); assert(index == PARTITION_ID_INITIAL); convert_geometry_objects(); // - update flags that tell whether we have reactions for all volume/surface species // and also update molecule type compartment flag // - must be done after geometry object conversions because // surface classes might have been defined world->get_all_rxns().update_all_mols_and_mol_type_compartments(); // we need to schedule the initial release for surfaces before the other releases world->create_initial_surface_region_release_event(); convert_release_events(); convert_mol_or_rxn_count_events_and_init_counting_flags(); convert_viz_output_events(); // beside of the event checking for check ctrl-c, it is a periodic event for each // iteration so, in the following call, where time maybe be moved time for scheduler // when a checkpoint is resumed, we always end up at the starting iteration and do not skip it // so that later we are not inserting events to the past (from the scheduler's point of view) add_ctrl_c_termination_event(); // sets up data loaded from checkpoint convert_initial_iteration_and_time_and_move_scheduler_time(); // some general checks if (world->config.rxn_radius_3d * POS_SQRT2 >= world->config.subpart_edge_length / 2) { throw ValueError(S("Reaction radius multiplied by sqrt(2) ") + to_string(world->config.rxn_radius_3d * world->config.length_unit * POS_SQRT2) + " must be less than half of subpartition edge length " + to_string(world->config.subpart_edge_length * world->config.length_unit / 2) + ". " + "Increase the model's " + NAME_CONFIG + "." + NAME_SUBPARTITION_DIMENSION + "."); } check_all_mol_types_have_diffusion_const(); } species_id_t MCell4Converter::get_species_id_for_complex( API::Complex& ci, const std::string error_msg, const bool check_orientation) { // check that the complex instance if fully qualified BNG::Cplx bng_ci = bng_converter.convert_complex(ci, true, !check_orientation); if (!bng_ci.is_fully_qualified()) { // TODO: add test throw ValueError( error_msg + ": " + NAME_COMPLEX + " '" + bng_ci.to_str() + "' must be fully qualified " + "(all components must be present and their state set)."); } // check that all used elementary molecule types have diffusion constant for (const BNG::ElemMol& em: bng_ci.elem_mols) { const BNG::ElemMolType& emt = world->bng_engine.get_data().get_elem_mol_type(em.elem_mol_type_id); if (emt.D == FLT_INVALID) { throw RuntimeError("Molecule type '" + emt.name + "' does not have its diffusion constant specified."); } } // we need to define species for our complex instance BNG::Species s = BNG::Species( bng_ci, world->bng_engine.get_data(), world->bng_engine.get_config() ); return world->bng_engine.get_all_species().find_or_add(s); } species_id_t MCell4Converter::get_species_id( API::Species& s, const std::string class_name, const std::string object_name) { if (s.species_id != SPECIES_ID_INVALID) { return s.species_id; } else { // we fist need to create a complex instance from our species API::Complex* s_as_cplx_inst = dynamic_cast(&s); return get_species_id_for_complex(*s_as_cplx_inst, class_name + " " + object_name); } } void MCell4Converter::get_geometry_bounding_box(Vec3& llf, Vec3& urb) { llf = Vec3(DBL_GIGANTIC); urb = Vec3(-DBL_GIGANTIC); for (std::shared_ptr& o: model->geometry_objects) { // go through all vertices for (auto& vert: o->vertex_list) { Vec3 v(vert); if (v.x < llf.x) { llf.x = v.x; } if (v.x > urb.x) { urb.x = v.x; } if (v.y < llf.y) { llf.y = v.y; } if (v.y > urb.y) { urb.y = v.y; } if (v.z < llf.z) { llf.z = v.z; } if (v.z > urb.z) { urb.z = v.z; } } } } void MCell4Converter::convert_simulation_setup() { // notifications and reports const API::Notifications& notifications = model->notifications; world->config.rxn_and_species_report = notifications.rxn_and_species_report; world->config.iteration_report = notifications.iteration_report; world->config.wall_overlap_report = notifications.wall_overlap_report; world->config.simulation_stats_every_n_iterations = notifications.simulation_stats_every_n_iterations; world->config.notifications.bng_verbosity_level = notifications.bng_verbosity_level; world->config.notifications.rxn_probability_changed = notifications.rxn_probability_changed; // warnings const API::Warnings& warnings = model->warnings; assert(warnings.high_reaction_probability != WarningLevel::ERROR); world->config.warnings.warn_on_bimol_rxn_probability_over_05_less_1 = warnings.high_reaction_probability == WarningLevel::WARNING; world->config.molecule_placement_failure = warnings.molecule_placement_failure; // config const API::Config& config = model->config; world->total_iterations = config.total_iterations; world->config.time_unit = config.time_step; world->config.use_bng_units = config.use_bng_units; world->config.initial_time = config.initial_time; world->config.initial_iteration = config.initial_iteration; pos_t grid_density = config.surface_grid_density; world->config.grid_density = grid_density; pos_t length_unit = 1/sqrt_f(config.surface_grid_density); world->config.length_unit = length_unit; if (is_set(config.interaction_radius)) { // NOTE: mcell3 does not convert the unit of the interaction radius in parser which // seems a bit weird world->config.rxn_radius_3d = config.interaction_radius / length_unit; } else { world->config.rxn_radius_3d = world->config.get_default_rxn_radius_3d(); } if (is_set(config.intermembrane_interaction_radius)) { // NOTE: mcell3 does not convert the unit of the interaction radius in parser which // seems a bit weird world->config.intermembrane_rxn_radius_3d = config.intermembrane_interaction_radius / length_unit; } else { world->config.intermembrane_rxn_radius_3d = (1.0 / sqrt_f(MY_PI * grid_density)) / length_unit; } pos_t vacancy_search_dist = config.vacancy_search_distance / length_unit; // Convert units world->config.vacancy_search_dist2 = vacancy_search_dist * vacancy_search_dist; // and take square world->config.randomize_smol_pos = !config.center_molecules_on_grid; world->config.sort_mols_by_subpart = config.sort_molecules; world->config.check_overlapped_walls = config.check_overlapped_walls; world->config.initial_seed = config.seed; rng_init(&world->rng, world->config.initial_seed); Vec3 llf, urb; get_geometry_bounding_box(llf, urb); Vec3 llf_w_margin = llf - Vec3(PARTITION_EDGE_EXTRA_MARGIN_UM); Vec3 urb_w_margin = urb + Vec3(PARTITION_EDGE_EXTRA_MARGIN_UM); pos_t llf_partition_dimension_diff = 0; // TODO: simplify this code to set origin and partition dimensions pos_t auto_partition_dimension = max3(urb_w_margin) - min3(llf_w_margin); bool auto_origin_set = false; if (!is_set(config.initial_partition_origin) && auto_partition_dimension > config.partition_dimension) { cout << "Info: Value of " << NAME_CLASS_MODEL << "." << NAME_CONFIG << "." << NAME_PARTITION_DIMENSION << " " << config.partition_dimension << " is smaller than the automatically determined value " << auto_partition_dimension << ", using the automatic value.\n"; world->config.partition_edge_length = auto_partition_dimension / length_unit; // we need to move the origin, if specified, by half of this increment llf_partition_dimension_diff = -(world->config.partition_edge_length - (config.partition_dimension / length_unit)) / 2; // place origin to the llf of the bounding box world->config.partition0_llf = (llf - Vec3(PARTITION_EDGE_EXTRA_MARGIN_UM)) / Vec3(length_unit); auto_origin_set = true; } else if (is_set(config.initial_partition_origin) && (config.initial_partition_origin[0] > llf_w_margin.x || config.initial_partition_origin[1] > llf_w_margin.y || config.initial_partition_origin[2] > llf_w_margin.z || config.initial_partition_origin[0] + config.partition_dimension < urb_w_margin.x || config.initial_partition_origin[1] + config.partition_dimension < urb_w_margin.y || config.initial_partition_origin[2] + config.partition_dimension < urb_w_margin.z ) ) { Vec3 origin(config.initial_partition_origin); Vec3 llf_diff(0); if (origin.x > llf_w_margin.x) { llf_diff.x = origin.x - llf_w_margin.x; } if (origin.y > llf_w_margin.y) { llf_diff.y = origin.y - llf_w_margin.y; } if (origin.z > llf_w_margin.z) { llf_diff.z = origin.z - llf_w_margin.z; } llf_partition_dimension_diff = max3(llf_diff); Vec3 opposite(Vec3(config.initial_partition_origin) + Vec3(config.partition_dimension)); Vec3 urb_diff(0); if (opposite.x < urb_w_margin.x) { urb_diff.x = urb_w_margin.x - opposite.x; } if (opposite.y < urb_w_margin.y) { urb_diff.y = urb_w_margin.y - opposite.y; } if (opposite.z < urb_w_margin.z) { urb_diff.z = urb_w_margin.z - opposite.z; } pos_t max_urb_diff = max3(urb_diff); world->config.partition_edge_length = (config.partition_dimension + llf_partition_dimension_diff + max_urb_diff)/ length_unit; // we need to move the origin, if specified, by half of this increment cout << "Info: Value of " << NAME_INITIAL_PARTITION_ORIGIN << " " << origin << " does not provide enough margin" " for model's geometry bounding box lower, left, front point " << llf << " " " and upper, right, back " << urb << "." " Moving the origin to " << origin - Vec3(llf_partition_dimension_diff) << " and increasing " << NAME_PARTITION_DIMENSION << " to " << world->config.partition_edge_length * length_unit << ".\n"; } else { world->config.partition_edge_length = config.partition_dimension / length_unit; } if (is_set(config.initial_partition_origin)) { // origin set manually world->config.partition0_llf = (Vec3(config.initial_partition_origin) - Vec3(llf_partition_dimension_diff)) / Vec3(length_unit); } else if (!auto_origin_set) { // place the partition to the center world->config.partition0_llf = -Vec3(world->config.partition_edge_length) / Vec3(2); } // align the origin to a multiple of subpartition length pos_t sp_len = config.subpartition_dimension / length_unit; uint tentative_subparts = world->config.partition_edge_length / sp_len; if (tentative_subparts > MAX_SUBPARTS_PER_PARTITION) { cout << "Info: Approximate number of subpartitions " << tentative_subparts << " is too high, lowering it to a limit of " << MAX_SUBPARTS_PER_PARTITION << ".\n"; sp_len = world->config.partition_edge_length / MAX_SUBPARTS_PER_PARTITION; } Vec3 orig_origin = world->config.partition0_llf; world->config.partition0_llf = floor_to_multiple_allow_negative_p(orig_origin, sp_len); // enlarge the partition by size we moved it in order to be aligned Vec3 llf_moved = orig_origin - world->config.partition0_llf; pos_t partition_edge_length_enlarged = world->config.partition_edge_length + max3(llf_moved); world->config.partition_edge_length = ceil_to_multiple_p(partition_edge_length_enlarged, sp_len); world->config.num_subparts_per_partition_edge = round_f(world->config.partition_edge_length / sp_len); // this option in MCell3 was removed in MCell4 world->config.use_expanded_list = true; // TODO: check that the values are higher than 1 world->config.rxn_class_cleanup_periodicity = config.reaction_class_cleanup_periodicity; world->config.species_cleanup_periodicity = config.species_cleanup_periodicity; world->config.molecules_order_random_shuffle_periodicity = config.molecules_order_random_shuffle_periodicity; world->config.memory_limit_gb = config.memory_limit_gb; world->config.continue_after_sigalrm = config.continue_after_sigalrm; // compute other constants and initialize reporting (if enabled) world->config.init(); } void MCell4Converter::convert_elementary_molecule_types() { BNG::BNGData& bng_data = world->bng_engine.get_data(); for (std::shared_ptr& api_mt: model->elementary_molecule_types) { bng_converter.convert_elementary_molecule_type(*api_mt); } } void MCell4Converter::convert_species() { for (std::shared_ptr& s: model->species) { BNG::Species new_species(world->bng_engine.get_data()); new_species.name = s->name; if (is_set(s->diffusion_constant_3d)) { new_species.D = s->diffusion_constant_3d; new_species.set_is_vol(); } else if (is_set(s->diffusion_constant_2d)) { new_species.D = s->diffusion_constant_2d; new_species.set_is_surf(); } else if (BNG::is_species_superclass(new_species.name)) { // these values are not really used, they are initialized for comparisons new_species.D = 0; new_species.space_step = 0; new_species.time_step = 0; } else { // this is a declaration of a possibly complex molecule, all used elementary molecules must be defined if (is_set(s->custom_time_step) || is_set(s->custom_space_step)) { throw ValueError("Species declaration " + new_species.name + " must not use " + NAME_CUSTOM_TIME_STEP + " nor " + NAME_CUSTOM_SPACE_STEP + " because it is derived from " + "its elementary molecule types." ); } BNG::Cplx bng_cplx(&world->bng_engine.get_data()); new_species.elem_mols = bng_converter.convert_complex(*s, false, false).elem_mols; new_species.finalize_species(world->config); if (!new_species.is_fully_qualified()) { throw ValueError("Species declaration " + new_species.name + " must use a fully qualified BNGL string for initialization."); } // check that all used elementary molecule types have their diffusion constant specified for (const BNG::ElemMol& em: new_species.elem_mols) { const BNG::ElemMolType& emt = world->bng_engine.get_data().get_elem_mol_type(em.elem_mol_type_id); if (emt.D == FLT_INVALID) { throw ValueError(S("Neither ") + NAME_DIFFUSION_CONSTANT_2D + " nor " + NAME_DIFFUSION_CONSTANT_3D + " was set for an elementary molecule " + emt.name + " used in " + NAME_CLASS_SPECIES + " " + s->to_bngl_str() + "."); } } // register and remember which species we created species_id_t new_species_id = world->get_all_species().find_or_add(new_species); s->species_id = new_species_id; // we completely converted this declaration continue; } new_species.set_flag(BNG::SPECIES_MOL_FLAG_TARGET_ONLY, s->target_only); // default is false // FIXME: the MolType below is created correctly only for simple species release_assert(s->elementary_molecules.size() <= 1 && "TODO: Complex species"); for (auto& mi: s->elementary_molecules) { release_assert(mi->components.empty()); } // find elementary molecule type for our species // must exist because it was added in Subsystem::unify_and_register_elementary_molecule_types BNG::elem_mol_type_id_t mol_type_id = world->bng_engine.get_data().find_elem_mol_type_id(s->name); release_assert(mol_type_id != BNG::ELEM_MOL_TYPE_ID_INVALID); BNG::ElemMol mol_inst; mol_inst.elem_mol_type_id = mol_type_id; new_species.elem_mols.push_back(mol_inst); new_species.finalize_species(world->config); species_id_t new_species_id = world->get_all_species().find_or_add(new_species); // remember which species we created s->species_id = new_species_id; // and also set superclass id for container if needed if (new_species.name == ALL_MOLECULES) { world->get_all_species().set_all_molecules_ids(new_species_id, mol_type_id); } else if (new_species.name == ALL_VOLUME_MOLECULES) { world->get_all_species().set_all_volume_molecules_ids(new_species_id, mol_type_id); } else if (new_species.name == ALL_SURFACE_MOLECULES) { world->get_all_species().set_all_surface_molecules_ids(new_species_id, mol_type_id); } } } void MCell4Converter::convert_surface_class_rxn( API::SurfaceProperty& sp, const BNG::Species& surface_reactant) { if (sp.type == SurfacePropertyType::REACTIVE) { // no need to add any reaction because reactions are defined manually return; } BNG::Cplx affected_pattern = bng_converter.convert_complex(*sp.affected_complex_pattern, false, true); BNG::RxnRule rxn(&world->bng_engine.get_data()); rxn.name = affected_pattern.to_str() + "+" + surface_reactant.name; switch (sp.type) { case API::SurfacePropertyType::ABSORPTIVE: if (affected_pattern.is_surf()) { // we are using special type for surf + absorptive surf class rxn.type = BNG::RxnType::AbsorbRegionBorder; } else { // the type is standard for vol + absorptive surf class rxn.type = BNG::RxnType::Standard; } break; case API::SurfacePropertyType::CONCENTRATION_CLAMP: rxn.type = BNG::RxnType::Standard; rxn.set_flag(BNG::RXN_FLAG_CREATED_FOR_CONCENTRATION_CLAMP); break; case API::SurfacePropertyType::FLUX_CLAMP: rxn.type = BNG::RxnType::Reflect; rxn.set_flag(BNG::RXN_FLAG_CREATED_FOR_FLUX_CLAMP); break; case API::SurfacePropertyType::REFLECTIVE: rxn.type = BNG::RxnType::Reflect; break; case API::SurfacePropertyType::TRANSPARENT: rxn.type = BNG::RxnType::Transparent; break; default: throw ValueError("Invalid SurfaceProperty type for " + surface_reactant.name + "."); } // all these reactions happen always rxn.base_rate_constant = DBL_GIGANTIC; // any compartment of the affected_pattern.set_compartment_id(BNG::COMPARTMENT_ID_NONE); // NONE is ANY in rxns orientation_t orient = convert_api_orientation(sp.affected_complex_pattern->orientation, true); rxn.append_reactant(affected_pattern); rxn.append_reactant(surface_reactant); // copies the input reactant // this is the default orientation used for surfaces in a reaction rxn.reactants[1].set_orientation(ORIENTATION_UP); // add reaction and remember mapping sp.rxn_rule_id = world->get_all_rxns().add_and_finalize(rxn); } void MCell4Converter::convert_surface_classes() { for (std::shared_ptr& sc: model->surface_classes) { // each surface class is represented by a species BNG::Species sc_species(world->bng_engine.get_data()); sc_species.name = sc->name; sc_species.set_is_reactive_surface(); // sets steps to 0 sc_species.space_step = 0; sc_species.time_step = 0; sc_species.D = FLT_INVALID; // diffusion constant has no meaning for surface classes // we must add a complex instance as the single molecule type in the new species // define a molecule type with no components BNG::ElemMolType mol_type; mol_type.name = sc_species.name; // name of the mol type is the same as for our species mol_type.set_is_reactive_surface(); BNG::elem_mol_type_id_t mol_type_id = world->bng_engine.get_data().find_or_add_elem_mol_type(mol_type); BNG::ElemMol mol_inst; mol_inst.elem_mol_type_id = mol_type_id; sc_species.elem_mols.push_back(mol_inst); // we do not care about compartments for surface classes sc_species.set_compartment_id(BNG::COMPARTMENT_ID_NONE); sc_species.finalize_species(world->config, false); species_id_t new_species_id = world->get_all_species().find_or_add(sc_species); // remember which species we created sc->species_id = new_species_id; // and we also need to add a reaction for each property if (sc->properties.empty()) { convert_surface_class_rxn(*dynamic_pointer_cast(sc), sc_species); } else { for (shared_ptr& sp: sc->properties) { convert_surface_class_rxn(*sp, sc_species); } } } } void MCell4Converter::check_intermembrane_surface_reaction(const BNG::RxnRule& rxn) { if (rxn.reactants.size() !=2 || rxn.products.size() != 2) { throw ValueError("Intermembrane reaction must have exactly 2 reactants and 2 products, error for " + rxn.to_str() + "."); } if (!rxn.reactants[0].is_surf() || !rxn.reactants[1].is_surf() || !rxn.products[0].is_surf() || !rxn.products[1].is_surf()) { throw ValueError("Intermembrane reaction's reactants and products must be all surface complexes, error for " + rxn.to_str() + "."); } // TODO: check compartments? } void MCell4Converter::convert_rxns() { BNG::BNGData& bng_data = world->bng_engine.get_data(); for (std::shared_ptr& r: model->reaction_rules) { bool is_reversible = is_set(r->rev_rate); BNG::RxnRule rxn(&bng_data); rxn.type = BNG::RxnType::Standard; if (is_set(r->fwd_rate)) { rxn.base_rate_constant = r->fwd_rate; } else { assert(!r->variable_rate.empty()); for (auto& time_and_rate: r->variable_rate) { assert(time_and_rate.size() == 2); BNG::RxnRateInfo info; info.time = time_and_rate[0] / world->config.time_unit; info.rate_constant = time_and_rate[1]; rxn.base_variable_rates.push_back(info); } rxn.update_variable_rxn_rate(0, nullptr); } for (std::shared_ptr& rinst: r->reactants) { // convert to BNG::ComplexInstance using existing or new BNG::molecule_id BNG::Cplx reactant = bng_converter.convert_complex(*rinst, false, true); rxn.append_reactant(reactant); } for (std::shared_ptr& pinst: r->products) { // convert to BNG::ComplexInstance using existing or new BNG::molecule_id BNG::Cplx product = bng_converter.convert_complex(*pinst, false, true); rxn.append_product(product); } // sets also flags for reactants and products rxn.finalize(); if (r->is_intermembrane_surface_reaction) { check_intermembrane_surface_reaction(rxn); rxn.set_is_intermembrane_surf_rxn(); // orientation is ANY in this case, we might need to update it later // TODO: add warning if user specified explicit orientation rxn.reactants[0].set_orientation(ORIENTATION_NONE); rxn.reactants[1].set_orientation(ORIENTATION_NONE); } string error_msg = BNG::process_compartments_and_set_orientations(bng_data, rxn); if (error_msg != "") { throw ValueError(error_msg); } if (is_set(r->name)) { rxn.name = r->name; } else { // must be called after conversion rxn.set_automatic_name(false); } // reverse reaction BNG::RxnRule rxn_rev(&bng_data); if (is_reversible) { rxn_rev.type = BNG::RxnType::Standard; rxn_rev.base_rate_constant = r->rev_rate; rxn_rev.reactants = rxn.products; rxn_rev.products = rxn.reactants; if (r->is_intermembrane_surface_reaction) { rxn.set_is_intermembrane_surf_rxn(); rxn.reactants[0].set_orientation(ORIENTATION_NONE); rxn.reactants[1].set_orientation(ORIENTATION_NONE); } rxn_rev.finalize(); string error_msg = BNG::process_compartments_and_set_orientations(bng_data, rxn_rev); if (error_msg != "") { throw ValueError(error_msg); } if (is_set(r->rev_name)) { rxn.name = r->rev_name; } else { rxn.set_automatic_name(false); } } // add reaction(s) and remember mapping r->fwd_rxn_rule_id = world->get_all_rxns().add_and_finalize(rxn); if (is_reversible) { r->rev_rxn_rule_id = world->get_all_rxns().add_and_finalize(rxn_rev); } // the ReactionRule object also need the world pointer r->world = world; } } MCell::wall_index_t MCell4Converter::convert_wall_and_add_to_geom_object( const API::GeometryObject& src_obj, const uint side, MCell::Partition& p, MCell::GeometryObject& dst_obj) { assert(src_obj.wall_list[side].size() == 3); // TODO LATER: there is really no reason to add walls in two steps, // can be simplified MCell::Wall& wall = p.add_uninitialized_wall(world->get_next_wall_id()); wall.object_id = dst_obj.id; wall.object_index = dst_obj.index; wall.side = side; for (uint i = 0; i < VERTICES_IN_TRIANGLE; i++) { wall.vertex_indices[i] = src_obj.vertex_indices[src_obj.wall_list[side][i]]; } wall.initialize_wall_constants(p); // add wall to subpartitions dst_obj.wall_indices.push_back(wall.index); return wall.index; } void MCell4Converter::convert_initial_surface_releases( const std::vector>& api_releases, std::vector& mcell_releases ) { for (auto api_rel: api_releases) { species_id_t species_id = get_species_id_for_complex(*api_rel->complex, NAME_CLASS_INITIAL_SURFACE_RELEASE); orientation_t orientation = convert_api_orientation(api_rel->complex->orientation); if (is_set(api_rel->number_to_release)) { if (api_rel->number_to_release > (double)UINT32_MAX) { throw ValueError(S("Value for ") + NAME_NUMBER_TO_RELEASE + " of a " + NAME_CLASS_INITIAL_SURFACE_RELEASE + " " + to_string(api_rel->number_to_release) + " is too high, the maximum allowed is " + to_string(UINT32_MAX) + "."); } mcell_releases.push_back( InitialSurfaceReleases(species_id, orientation, true, (uint)api_rel->number_to_release) ); } else if (is_set(api_rel->density)) { mcell_releases.push_back( InitialSurfaceReleases(species_id, orientation, false, (double)api_rel->density) ); } else { assert(false); } } } void MCell4Converter::convert_concentration_clamp_release( const partition_id_t partition_id, const API::SurfaceClass& surface_class, const MCell::Region& mcell_region) { release_assert(surface_class.properties.empty() && "TODO"); ClampReleaseEvent* clamp_event = new ClampReleaseEvent(world); if (surface_class.type == SurfacePropertyType::CONCENTRATION_CLAMP) { clamp_event->type = ClampType::CONCENTRATION; } else if (surface_class.type == SurfacePropertyType::FLUX_CLAMP) { clamp_event->type = ClampType::FLUX; } else { assert(false); } // run each timestep clamp_event->event_time = 0; clamp_event->periodicity_interval = 1; // which species to clamp clamp_event->species_id = get_species_id_for_complex( *surface_class.affected_complex_pattern, S(NAME_CLASS_SURFACE_CLASS) + ", attribute " + NAME_AFFECTED_COMPLEX_PATTERN, false); // on which side clamp_event->orientation = convert_api_orientation(surface_class.affected_complex_pattern->orientation, true, true); assert(world->bng_engine.get_all_species().get(surface_class.species_id).is_reactive_surface()); clamp_event->surf_class_species_id = surface_class.species_id; clamp_event->concentration = surface_class.concentration; // walls where to release assert(!mcell_region.walls_and_edges.empty() && "Must be initialized"); // we are inserting the walls ordered by their wall index for (const auto& pair_wi_edges: mcell_region.walls_and_edges) { clamp_event->cumm_area_and_pwall_index_pairs.push_back( CummAreaPWallIndexPair(0, PartitionWallIndexPair(partition_id, pair_wi_edges.first))); } clamp_event->update_cumm_areas_and_scaling(); // and schedule world->scheduler.schedule_event(clamp_event); } MCell::region_index_t MCell4Converter::convert_surface_region( MCell::Partition& p, API::SurfaceRegion& surface_region, API::GeometryObject& o, MCell::GeometryObject& obj) { MCell::Region reg; reg.name = obj.name + "," + surface_region.name; reg.geometry_object_id = obj.id; // simply add all walls for (const int wall_in_object: surface_region.wall_indices) { wall_index_t wi = o.wall_indices[wall_in_object]; reg.add_wall_to_walls_and_edges(wi, false); } if (is_set(surface_region.surface_class)) { assert(surface_region.surface_class->species_id != SPECIES_ID_INVALID); reg.species_id = surface_region.surface_class->species_id; // define releases for concentration clamp if (surface_region.surface_class->is_clamp()) { convert_concentration_clamp_release(p.id, *surface_region.surface_class, reg); } } if (is_set(surface_region.initial_surface_releases)) { convert_initial_surface_releases( surface_region.initial_surface_releases, reg.initial_region_molecules ); } reg.init_surface_region_edges(p); reg.id = world->get_next_region_id(); surface_region.region_id = reg.id; region_index_t ri = p.add_region_and_set_its_index(reg); return ri; } void MCell4Converter::convert_geometry_objects() { MCell::Partition& p = world->get_partition(PARTITION_ID_INITIAL); // only partition 0 is supported for now for (std::shared_ptr& o: model->geometry_objects) { // set surface region parents for (auto& sr: o->surface_regions) { // not using shared pointers here, any attempt so far resulted in bad_weak_ptr exception // this is safe because the geometry object (parent) has a reference to the surface region sr->parent = o.get(); } MCell::GeometryObject& obj = p.add_uninitialized_geometry_object(world->get_next_geometry_object_id()); obj.name = o->name; obj.parent_name = ""; // empty, we do not really care o->geometry_object_id = obj.id; if (is_set(o->initial_color)) { // set default color obj.default_color = o->initial_color->rgba; } // vertices // remember the "offset" from the first vertex in the target partition o->first_vertex_index = p.get_geometry_vertex_count(); for (auto& v: o->vertex_list) { // add to partition and remember its index // must use rcp_length_unit to be identical to mcell4 with mdl vertex_index_t vi = p.add_geometry_vertex(Vec3(v[0], v[1], v[2]) * Vec3(world->config.rcp_length_unit)); o->vertex_indices.push_back(vi); } // walls (validity of indices is checked in API::GeometryObject::check_semantics) o->first_wall_index = p.get_walls().size(); for (size_t i = 0; i < o->wall_list.size(); i++) { wall_index_t wi = convert_wall_and_add_to_geom_object(*o, i, p, obj); o->wall_indices.push_back(wi); } // initialize edges obj.initialize_neighboring_walls_and_their_edges(p); vector region_indices; // region "ALL" MCell::Region reg_all; reg_all.init_from_whole_geom_object(obj); reg_all.id = world->get_next_region_id(); // TODO: move surf class handling code to a function and share with convert_surface_region if (is_set(o->surface_class)) { if (o->surface_class->species_id == SPECIES_ID_INVALID) { throw RuntimeError("Geometry object " + o->name + " is using surface class " + o->surface_class->name + " but this surface class was not added to the model. To add it, use method Model.add_surface_class or Subsystem.add_surface_class."); } reg_all.species_id = o->surface_class->species_id; // define releases for concentration clamp if (o->surface_class->is_clamp()) { convert_concentration_clamp_release(p.id, *o->surface_class, reg_all); } } if (is_set(o->initial_surface_releases)) { convert_initial_surface_releases(o->initial_surface_releases, reg_all.initial_region_molecules); } region_index_t ri_all = p.add_region_and_set_its_index(reg_all); region_indices.push_back(ri_all); // we must remember that this region represents the whole object obj.encompassing_region_id = reg_all.id; // regions from surface areas // mcell3 stores the regions in reverse, so we can too... // (maybe change this in the mcell3 converter) for (int k = (int)o->surface_regions.size() - 1; k >= 0; k--) { std::shared_ptr surface_region = o->surface_regions[k]; region_index_t ri = convert_surface_region(p, *surface_region, *o, obj); region_indices.push_back(ri); // also set colors if they were specified if (is_set(surface_region->initial_color)) { rgba_t color = surface_region->initial_color->rgba; for (wall_index_t wi: surface_region->wall_indices) { obj.wall_specific_colors[o->get_partition_wall_index(wi)] = color; } } } // also set regions for walls for (wall_index_t wi: obj.wall_indices) { MCell::Wall& w = p.get_wall(wi); for (region_index_t ri: region_indices) { MCell::Region& reg = p.get_region(ri); if (reg.walls_and_edges.count(wi) == 1) { w.regions.insert(ri); obj.regions.insert(ri); } } } // set MCell:GeometryObject compartment ids if (o->is_bngl_compartment) { assert(o->vol_compartment_id != BNG::COMPARTMENT_ID_INVALID && o->vol_compartment_id != BNG::COMPARTMENT_ID_NONE ); obj.vol_compartment_id = o->vol_compartment_id; if (is_set(o->surface_compartment_name)) { assert(o->surf_compartment_id != BNG::COMPARTMENT_ID_INVALID); obj.surf_compartment_id = o->surf_compartment_id; } } } for (MCell::GeometryObject& obj: p.get_geometry_objects()) { obj.initialize_is_fully_transparent(p); } // check overlapped walls // uses random generator state if (world->config.check_overlapped_walls) { bool ok = world->check_for_overlapped_walls(); if (!ok) { throw ValueError("Walls in geometry overlap, more details were printed in the previous message."); } } // check that no overlapped object is a BNGL compartment for (MCell::GeometryObject& obj: p.get_geometry_objects()) { if (obj.has_overlapped_walls && obj.vol_compartment_id != BNG::COMPARTMENT_ID_NONE) { throw ValueError("Overlapped geometry object " + obj.name + " is set to be a BNGL compartment, this is not allowed."); } } p.finalize_walls(); } void MCell4Converter::check_surface_compartment_name_collision(const std::string& surface_compartment_name) { for (std::shared_ptr& o: model->geometry_objects) { for (std::shared_ptr& s: o->surface_regions) { // o->wall_indices contains all the walls (with values counted from 0), so if size is the same // the contents are the same if (s->name == surface_compartment_name && o->wall_list.size() != s->wall_indices.size()) { throw RuntimeError("Geometry object's " + o->name + " surface region " + s->name + " uses the same name " "as a compartment, but the surface region does not represent the whole surface of the object. " "This is not allowed because it would lead to inconsistencies when comparing to BNGL variant of the model. " "Please use a different name either for the surface region or for the compartment."); } } } } void MCell4Converter::convert_compartments() { // we determine the hierarchy of compartments here since all // we have are geometry objects // volume of compartments is not set BNG::BNGData& bng_data = world->bng_engine.get_data(); vector> compartment_objects; for (std::shared_ptr& o: model->geometry_objects) { if (o->is_bngl_compartment) { compartment_objects.push_back(o); } } // set hierarchy of compartments set_parent_and_children_compartments(compartment_objects); for (std::shared_ptr& o: compartment_objects) { BNG::Compartment bng_comp3d; bng_comp3d.name = o->name; bng_comp3d.is_3d = true; BNG::compartment_id_t comp3d_id = bng_data.add_compartment(bng_comp3d); o->vol_compartment_id = comp3d_id; // unlike as in BNG, we do not require that the only child of a 3d compartment is 2d compartment, // 2d compartments can be skipped completely if (is_set(o->surface_compartment_name)) { // check that there is no SurfaceRegion with this name // this might be supported one day but now the surface compartment must be the whole // surface of the object, reported in redmine as #28 check_surface_compartment_name_collision(o->surface_compartment_name); BNG::Compartment bng_comp2d; bng_comp2d.name = o->surface_compartment_name; bng_comp2d.is_3d = false; // the only child of a 2D compartment is the 3D compartment it encompasses bng_comp2d.children_compartments.insert(comp3d_id); BNG::compartment_id_t comp2d_id = bng_data.add_compartment(bng_comp2d); o->surf_compartment_id = comp2d_id; // if a 2d compartment is defined, it is the parent of the 3D compartment bng_data.get_compartment(comp3d_id).parent_compartment_id = comp2d_id; } } // now define their parents and children for (std::shared_ptr& o: compartment_objects) { BNG::Compartment* bng_comp3d; bng_comp3d = bng_data.find_compartment(o->name); release_assert(bng_comp3d != nullptr); for (const std::shared_ptr& child: o->child_compartments) { BNG::Compartment* bng_comp_child; // the direct child is the 2d compartment if defined, 3d compartment otherwise if (is_set(child->surface_compartment_name)) { bng_comp_child = bng_data.find_compartment(child->surface_compartment_name); release_assert(bng_comp_child != nullptr); } else { bng_comp_child = bng_data.find_compartment(child->name); release_assert(bng_comp_child != nullptr); } bng_comp3d->children_compartments.insert(bng_comp_child->id); bng_comp_child->parent_compartment_id = bng_comp3d->id; } } } static MCell::RegionExprOperator convert_region_node_type(API::RegionNodeType t) { switch (t) { case API::RegionNodeType::LEAF_GEOMETRY_OBJECT: return MCell::RegionExprOperator::LEAF_GEOMETRY_OBJECT; case API::RegionNodeType::LEAF_SURFACE_REGION: return MCell::RegionExprOperator::LEAF_SURFACE_REGION; case API::RegionNodeType::UNION: return MCell::RegionExprOperator::UNION; case API::RegionNodeType::DIFFERENCE: return MCell::RegionExprOperator::DIFFERENCE; case API::RegionNodeType::INTERSECT: return MCell::RegionExprOperator::INTERSECT; default: assert(false); return MCell::RegionExprOperator::INVALID; } } RegionExprNode* MCell4Converter::convert_region_expr_recursively( const shared_ptr& region, const bool is_vol, const bool release_not_count, // for error message purposes MCell::RegionExpr& region_expr ) { string msg_detail; if (release_not_count) { msg_detail = S(" referenced by a ") + NAME_CLASS_RELEASE_SITE + " object"; } else { msg_detail = S(" referenced by a ") + NAME_CLASS_COUNT + " object"; } assert(is_set(region)); if (region->node_type == RegionNodeType::LEAF_GEOMETRY_OBJECT) { shared_ptr geometry_object = dynamic_pointer_cast(region); assert(is_set(geometry_object)); if (is_vol) { if (geometry_object->geometry_object_id == GEOMETRY_OBJECT_ID_INVALID) { throw RuntimeError("Could not find geometry object with name " + geometry_object->name + msg_detail + ", it might not have been added to the model."); } return region_expr.create_new_expr_node_leaf_geometry_object(geometry_object->geometry_object_id); } else { const MCell::Region* reg = world->find_region_by_name(geometry_object->name + MCell::REGION_ALL_SUFFIX_W_COMMA); if (reg == nullptr) { throw RuntimeError("Could not find region representing object with name " + geometry_object->name + msg_detail + ", it might not have been added to the model."); } return region_expr.create_new_expr_node_leaf_surface_region(reg->id); } } else if (region->node_type == RegionNodeType::LEAF_SURFACE_REGION) { if (is_vol) { // TODO: message needs more details string msg; if (release_not_count) { msg = "Cannot release volume molecule onto a surface region ."; } else { msg = "Cannot count volume molecules or reactions on a surface region ."; } throw RuntimeError(msg + "The problematic region is:\n" + region->to_str()); } shared_ptr surface_region = dynamic_pointer_cast(region); assert(is_set(surface_region)); const MCell::Region* reg = world->find_region_by_name(surface_region->parent->name + "," + surface_region->name); if (reg == nullptr) { throw RuntimeError("Could not find region with name " + surface_region->parent->name + msg_detail + ", it might not have been added to the model."); } return region_expr.create_new_expr_node_leaf_surface_region(reg->id); } else { return region_expr.create_new_region_expr_node_op( convert_region_node_type(region->node_type), convert_region_expr_recursively(region->left_node, is_vol, release_not_count, region_expr), convert_region_expr_recursively(region->right_node, is_vol, release_not_count, region_expr) ); } } void MCell4Converter::convert_rel_site_region_expr(API::ReleaseSite& rel_site, MCell::ReleaseEvent* rel_event) { const BNG::Species& species = world->get_all_species().get(rel_event->species_id); if (rel_site.shape == Shape::COMPARTMENT && !is_set(rel_site.complex->compartment_name)) { // this should not happen throw RuntimeError("Compartment for release site " + rel_site.name + " was not set."); } if (rel_site.shape == Shape::REGION_EXPR) { if (!is_set(rel_site.region)) { throw RuntimeError("Region for release site " + rel_site.name + " was not set."); } rel_event->region_expr.root = convert_region_expr_recursively( rel_site.region, species.is_vol(), true, rel_event->region_expr); // add intersection with compartment if this is a volume compartment if (is_set(rel_site.complex->compartment_name)) { // TODO: how to handle surface compartment intersection? // obj is nullptr when this is not a volume compartment const auto& obj = model->find_volume_compartment_object(rel_site.complex->compartment_name); if (is_set(obj)) { auto region = model->get_compartment_region(rel_site.complex->compartment_name); RegionExprNode* compartment_region = convert_region_expr_recursively( region, species.is_vol(), true, rel_event->region_expr); // overwrite region expr with intersection with compartment rel_event->region_expr.root = rel_event->region_expr.create_new_region_expr_node_op( RegionExprOperator::INTERSECT, rel_event->region_expr.root, compartment_region ); } } } else if (rel_site.shape == Shape::COMPARTMENT) { if (!is_set(rel_site.complex->compartment_name)) { throw RuntimeError("Compartment for release site " + rel_site.name + " was not set."); } // make region that represents the compartment auto region = model->get_compartment_region(rel_site.complex->compartment_name); if (!is_set(region)) { throw RuntimeError("Compartment " + rel_site.complex->compartment_name + " for " + rel_site.name + " was not found."); } rel_event->region_expr.root = convert_region_expr_recursively( region, species.is_vol(), true, rel_event->region_expr); } // also set llf and urb // TODO: this does not check anything yet bool ok = Geometry::compute_region_expr_bounding_box(world, rel_event->region_expr.root, rel_event->region_llf, rel_event->region_urb); if (!ok) { throw RuntimeError("Region for releases specified by " + rel_event->region_expr.root->to_string(world) + " is not closed."); } } void MCell4Converter::convert_molecule_list( const std::vector>& molecule_list, const std::string& rel_site_name, MCell::ReleaseEvent* rel_event) { for (auto& item: molecule_list) { MCell::SingleMoleculeReleaseInfo info; info.species_id = get_species_id_for_complex( *item->complex, S(NAME_CLASS_RELEASE_SITE) + " '" + rel_site_name + "'"); assert(item->location.size() == 3); info.pos.x = item->location[0] * world->config.rcp_length_unit; info.pos.y = item->location[1] * world->config.rcp_length_unit; info.pos.z = item->location[2] * world->config.rcp_length_unit; bool is_vol = world->get_all_species().get(info.species_id).is_vol(); info.orientation = convert_api_orientation(item->complex->orientation, true, is_vol); // not set is not allowed if (world->get_all_species().get(info.species_id).is_vol() && rel_event->orientation != ORIENTATION_NONE && rel_event->orientation != ORIENTATION_NOT_SET) { throw ValueError( S(NAME_CLASS_RELEASE_SITE) + " " + rel_site_name + " releases a volume molecule but orientation is set."); } rel_event->molecule_list.push_back(info); } } MCell::ReleaseEvent* MCell4Converter::convert_single_release_event( const std::shared_ptr& r) { MCell::ReleaseEvent* rel_event = new ReleaseEvent(world); rel_event->release_site_name = r->name; if (!is_set(r->molecule_list)) { assert(is_set(r->complex)); rel_event->species_id = get_species_id_for_complex( *r->complex, S(NAME_CLASS_RELEASE_SITE) + " '" + r->name + "'"); bool is_vol = world->get_all_species().get(rel_event->species_id).is_vol(); rel_event->orientation = convert_api_orientation(r->complex->orientation, true, is_vol); if (world->get_all_species().get(rel_event->species_id).is_surf() && rel_event->orientation != ORIENTATION_UP && rel_event->orientation != ORIENTATION_DOWN && rel_event->orientation != ORIENTATION_NONE) { // none = random orientation throw ValueError( S(NAME_CLASS_RELEASE_SITE) + " " + r->name + " releases a surface molecule but orientation is not set to a valid value."); } if (is_vol && rel_event->orientation != ORIENTATION_NONE && rel_event->orientation != ORIENTATION_NOT_SET) { throw ValueError( S(NAME_CLASS_RELEASE_SITE) + " " + r->name + " releases a volume molecule but orientation is set."); } } rel_event->delay = r->release_time / world->config.time_unit; // release pattern if (is_set(r->release_pattern)) { API::ReleasePattern rp = *r->release_pattern; rel_event->number_of_trains = rp.number_of_trains; rel_event->train_interval = rp.train_interval / world->config.time_unit; rel_event->train_duration = rp.train_duration / world->config.time_unit; rel_event->release_interval = rp.release_interval / world->config.time_unit; } // release_number_method if (is_set(r->number_to_release)) { if (r->number_to_release > (double)UINT32_MAX) { throw ValueError(S("Value ") + to_string(r->number_to_release) + " for " + NAME_NUMBER_TO_RELEASE + " of a " + NAME_CLASS_RELEASE_SITE + " '" + r->name + "' is too high, the maximum allowed is " + to_string(UINT32_MAX) + "."); } rel_event->release_number_method = ReleaseNumberMethod::CONST_NUM; rel_event->release_number = r->number_to_release; } else if (is_set(r->density)) { rel_event->release_number_method = ReleaseNumberMethod::DENSITY_NUM; rel_event->concentration = r->density; } else if (is_set(r->concentration)) { rel_event->release_number_method = ReleaseNumberMethod::CONCENTRATION_NUM; rel_event->concentration = r->concentration; } else if (is_set(r->molecule_list)) { rel_event->release_number_method = ReleaseNumberMethod::CONST_NUM; convert_molecule_list(r->molecule_list, r->name, rel_event); } else { throw RuntimeError( S("The only supported release number type now is constant number specified with ") + NAME_NUMBER_TO_RELEASE + "." ); } // release_shape switch (r->shape) { case Shape::SPHERICAL: rel_event->release_shape = ReleaseShape::SPHERICAL; assert(r->location.size() == 3); rel_event->location = Vec3(r->location) * world->config.rcp_length_unit; rel_event->diameter = r->site_diameter * world->config.rcp_length_unit; break; case Shape::REGION_EXPR: case Shape::COMPARTMENT: { rel_event->release_shape = ReleaseShape::REGION; convert_rel_site_region_expr(*r, rel_event); bool ok = rel_event->initialize_walls_for_release(); if (!ok) { throw RuntimeError("Only simple surface regions are supported for surface releases currently, error for " + r->name + "."); } } break; case Shape::LIST: rel_event->release_shape = ReleaseShape::LIST; rel_event->diameter = r->site_diameter * world->config.rcp_length_unit; break; rel_event->release_shape = ReleaseShape::REGION; break; default: // should be caught earlier throw RuntimeError(S("The only supported shapes now are ") + NAME_EV_SPHERICAL + ", " + NAME_EV_REGION_EXPR + ", " + NAME_EV_COMPARTMENT + " and " + NAME_EV_LIST + "."); } rel_event->release_probability = r->release_probability; return rel_event; } void MCell4Converter::convert_release_events() { for (std::shared_ptr& r: model->release_sites) { MCell::ReleaseEvent* rel_event = convert_single_release_event(r); // schedule it rel_event->update_event_time_for_next_scheduled_time(); world->scheduler.schedule_event(rel_event); } } static void set_geometry_objects_are_used_in_mol_rxn_counts_recursively( World* world, MCell::RegionExprNode* node) { if (node->op == RegionExprOperator::LEAF_GEOMETRY_OBJECT) { world->get_geometry_object(node->geometry_object_id).set_is_used_in_mol_rxn_counts(); } else if (node->has_binary_op()) { set_geometry_objects_are_used_in_mol_rxn_counts_recursively(world, node->left); set_geometry_objects_are_used_in_mol_rxn_counts_recursively(world, node->right); } else { release_assert(false && "Surface regions cannot be used here"); } } // appends new term to the vector terms, does not clear it void MCell4Converter::convert_count_term_leaf_and_init_counting_flags( const std::shared_ptr ct, const int sign, MolOrRxnCountTermVector& terms ) { MCell::MolOrRxnCountTerm res; res.sign_in_expression = sign; assert(is_set(ct)); assert(ct->node_type == API::ExprNodeType::LEAF); if (is_set(ct->species_pattern) || is_set(ct->molecules_pattern)) { if (is_set(ct->species_pattern)) { res.species_pattern_type = SpeciesPatternType::SpeciesPattern; res.species_molecules_pattern = bng_converter.convert_complex(*ct->species_pattern, true); } else { res.species_pattern_type = SpeciesPatternType::MoleculesPattern; res.species_molecules_pattern = bng_converter.convert_complex(*ct->molecules_pattern, true); } // to maintain compatibility with BioNetGen counting, we must remove the primary compartment // from elementary molecules and keep it separately, e.g.: pattern @PM:V(s!1).S(v!1) // which is after being read from BNGL in reality V(s!1)@PM.S(v!1)@PM must match // V(s!1)@CP.S(v!1)@PM, therefore if we make the pattern to be V(s!1).S(v!1) + PM, it will match res.primary_compartment_id = res.species_molecules_pattern.get_primary_compartment_id(); res.species_molecules_pattern.remove_compartment_from_elem_mols(res.primary_compartment_id); bool is_vol = res.species_molecules_pattern.is_vol(); string name = res.species_molecules_pattern.to_str(); res.orientation = res.species_molecules_pattern.get_orientation(); if (is_set(ct->region)) { if (is_vol) { res.type = MCell::CountType::EnclosedInVolumeRegion; res.region_expr.root = convert_region_expr_recursively(ct->region, is_vol, false, res.region_expr); // and also mark the objects that we are counting molecules inside set_geometry_objects_are_used_in_mol_rxn_counts_recursively(world, res.region_expr.root); } else { // surf mols res.type = MCell::CountType::PresentOnSurfaceRegion; res.region_expr.root = convert_region_expr_recursively(ct->region, is_vol, false, res.region_expr); // these are only surface regions and there is no need to set that they are counted } } else { res.type = MCell::CountType::EnclosedInWorld; } } else if (is_set(ct->reaction_rule)) { assert(!is_set(ct->reaction_rule->rev_rate)); res.rxn_rule_id = ct->reaction_rule->fwd_rxn_rule_id; // is this a surface rxn? -> at least one of the reactants is a surface mol BNG::RxnRule* rxn = world->get_all_rxns().get(res.rxn_rule_id); if (is_set(ct->region)) { // TODO: what about counting reactions with surface classes? bool is_vol = !rxn->is_surf_rxn(); if (is_vol) { // volume reaction res.type = MCell::CountType::RxnCountInVolumeRegion; rxn->set_is_counted_in_volume_regions(); res.region_expr.root = convert_region_expr_recursively(ct->region, is_vol, false, res.region_expr); set_geometry_objects_are_used_in_mol_rxn_counts_recursively(world, res.region_expr.root); } else { // surface reaction res.type = MCell::CountType::RxnCountOnSurfaceRegion; rxn->set_is_counted_on_surface_regions(); res.region_expr.root = convert_region_expr_recursively(ct->region, is_vol, false, res.region_expr); // these are only surface regions and there is no need to set that they are counted } } else { res.type = MCell::CountType::RxnCountInWorld; rxn->set_is_counted_in_world(); } // remember the initial reactions count when resuming from a checkpoint, default is 0 res.initial_reactions_count = ct->initial_reactions_count; } else { assert(false); } terms.push_back(res); } void MCell4Converter::convert_count_terms_recursively( const std::shared_ptr count, // only for warning printouts const std::shared_ptr ct, const int sign, MCell::MolOrRxnCountItem& info ) { assert(is_set(ct)); if (ct->node_type == API::ExprNodeType::LEAF) { convert_count_term_leaf_and_init_counting_flags(ct, sign, info.terms); } else if (ct->node_type == API::ExprNodeType::ADD || ct->node_type == API::ExprNodeType::SUB) { if (is_set(ct->region)) { warns() << "Object " << NAME_CLASS_COUNT << " with " << NAME_NAME << " '" << count->name << "' and " << NAME_FILE_NAME << " '" << count->file_name << "' uses a " << NAME_CLASS_COUNT_TERM << " that has a " << NAME_REGION << " set. This region will be ignored because regions are allowed" << " only for leaf " << NAME_CLASS_COUNT_TERM << " objects.\n" << "This is the problematic " << NAME_CLASS_COUNT << " object:\n" << count->to_str() << "\n" << "-- end of warning message reporting ignored region --\n"; } convert_count_terms_recursively(count, ct->left_node, sign, info); int next_sign; if (ct->node_type == API::ExprNodeType::SUB) { next_sign = -sign; } else { next_sign = sign; } convert_count_terms_recursively(count, ct->right_node, next_sign, info); } else { // cannot really happen throw RuntimeError("Invalid node_type in CountTerm."); } } void MCell4Converter::convert_mol_or_rxn_count_events_and_init_counting_flags() { // collect counts with gdat format outputting to the same file std::map> gdat_filename_to_count_indices; for (uint i = 0 ; i < model->counts.size(); i++) { std::shared_ptr& c = model->counts[i]; // set default file_name if not set if (!is_set(c->file_name)) { c->file_name = "./react_data/" + get_seed_dir_name(model->config.seed) + "/" + c->name + ".dat"; } if (c->output_format == CountOutputFormat::AUTOMATIC_FROM_EXTENSION) { throw RuntimeError(S(NAME_CLASS_COUNT) + "'s " + NAME_OUTPUT_FORMAT + " must not be " + NAME_ENUM_COUNT_OUTPUT_FORMAT + "." + NAME_EV_AUTOMATIC_FROM_EXTENSION + " when the model is initialized. " "The automatic detection should have already happened."); } if (c->output_format != CountOutputFormat::GDAT) { continue; } auto it = gdat_filename_to_count_indices.find(c->file_name); if (it == gdat_filename_to_count_indices.end()) { gdat_filename_to_count_indices[c->file_name].push_back(i); } else { it->second.push_back(i); } } // create GDAT output buffers std::map> gdat_count_name_to_buffer_id_and_column_index; for (const auto& pair_fname_indices: gdat_filename_to_count_indices) { // prepare names for this single gdat buffer // also check that the sampling interval is the same assert(pair_fname_indices.second.size() >= 1); double every_n_timesteps = model->counts[pair_fname_indices.second[0]]->every_n_timesteps; vector names; for (uint i: pair_fname_indices.second){ std::shared_ptr& c = model->counts[i]; if (c->every_n_timesteps != every_n_timesteps) { throw RuntimeError(S("When multiple ") + NAME_CLASS_COUNT + " objects output to the same .gdat file, " + "their sampling intervals set by " + NAME_EVERY_N_TIMESTEPS + " must be identical, error for " + c->name + "."); } names.push_back(c->name); } count_buffer_id_t buffer_id = world->create_gdat_count_buffer( pair_fname_indices.first, names, API::DEFAULT_COUNT_BUFFER_SIZE, model->config.append_to_count_output_data); // these are local column indices and must start from 0 for (uint i = 0; i < pair_fname_indices.second.size(); i++){ // however tpo get the count we must use the global index std::shared_ptr& c = model->counts[pair_fname_indices.second[i]]; gdat_count_name_to_buffer_id_and_column_index[c->name] = make_pair(buffer_id, i); } } // and convert the count objects for (const std::shared_ptr& c: model->counts) { MCell::MolOrRxnCountEvent* count_event = new MCell::MolOrRxnCountEvent(world); count_event->event_time = 0; count_event->periodicity_interval = round_f(c->every_n_timesteps + EPS); // create buffer vector column_names = {c->name}; count_buffer_id_t buffer_id; uint column_index; if (c->output_format == CountOutputFormat::DAT) { buffer_id = world->create_dat_count_buffer( c->file_name, API::DEFAULT_COUNT_BUFFER_SIZE, model->config.append_to_count_output_data); column_index = 0; } else { auto it = gdat_count_name_to_buffer_id_and_column_index.find(c->name); assert(it != gdat_count_name_to_buffer_id_and_column_index.end()); buffer_id = it->second.first; column_index = it->second.second; } MCell::MolOrRxnCountItem info(buffer_id, column_index); // process count terms convert_count_terms_recursively(c, c->expression, +1, info); info.multiplier = c->multiplier; // having multiple MolOrRxnCountInfo per MolOrRxnCountEvent // was useful for MCell3 conversion, however for pymcell4 each count is a separate event count_event->add_mol_count_item(info); // remember for get_current_value calls c->count_event = count_event; if (c->every_n_timesteps > 0) { // 0 means that the event won't be ever rescheduled, // if it were added to the schedule it would be executed once and then deleted world->scheduler.schedule_event(count_event); } else { world->add_unscheduled_count_event(count_event); } } } void MCell4Converter::convert_viz_output_events() { for (std::shared_ptr& v: model->viz_outputs) { MCell::VizOutputEvent* viz_event = new VizOutputEvent(world); viz_event->event_time = 0.0; viz_event->periodicity_interval = round_f(v->every_n_timesteps + EPS); viz_event->viz_mode = convert_viz_mode(v->mode); if (!is_set(v->output_files_prefix)) { v->output_files_prefix = "./viz_data/" + get_seed_dir_name(model->config.seed) + "/Scene"; } viz_event->file_prefix_name = v->output_files_prefix; if (is_set(v->species_list)) { for (std::shared_ptr& s: v->species_list) { viz_event->species_ids_to_visualize.insert( get_species_id(*s, NAME_CLASS_VIZ_OUTPUT, NAME_SPECIES_LIST)); } } else { viz_event->visualize_all_species = true; } world->scheduler.schedule_event(viz_event); } } // sets up data loaded from checkpoint, must be run after all events were added to the scheduler // rng state is set after model initialization in Model::initialize void MCell4Converter::convert_initial_iteration_and_time_and_move_scheduler_time() { if ((model->config.initial_iteration != 0) != (model->config.initial_time != 0)) { throw RuntimeError(S("Both ") + NAME_INITIAL_ITERATION + " and " + NAME_INITIAL_TIME + " must be either 0 or both must be greater than 0."); } world->stats.set_current_iteration(model->config.initial_iteration); // skips all events, fails if some periodic events are scheduled world->scheduler.skip_events_up_to_time( world->config.get_simulation_start_time() ); } void MCell4Converter::add_ctrl_c_termination_event() { MCell::CustomFunctionCallEvent* event = new CustomFunctionCallEvent(check_ctrl_c, world); event->event_time = world->config.get_simulation_start_time(); event->periodicity_interval = 1; world->scheduler.schedule_event(event); } void MCell4Converter::check_all_mol_types_have_diffusion_const() { for (const BNG::ElemMolType& mt: world->bng_engine.get_data().get_elem_mol_types()) { if (!mt.is_reactive_surface() && mt.D == FLT_INVALID) { throw RuntimeError("Molecule type '" + mt.name + "' does not have its diffusion constant specified."); } } } // must be called after world initialization void MCell4Converter::convert_after_init() { convert_rng_state(); convert_checkpointed_molecules(); } void MCell4Converter::convert_rng_state() { if (!is_set(model->config.initial_rng_state)) { return; // nothing to do } std::shared_ptr& src = model->config.initial_rng_state; rng_state& dst = world->rng; dst.randcnt = src->randcnt; dst.aa = src->aa; dst.bb = src->bb; dst.cc = src->cc; assert(RANDSIZ == RNG_SIZE); assert(src->randslr.size() == RNG_SIZE); std::copy(src->randslr.begin(), src->randslr.end(), dst.randrsl); assert(src->mm.size() == RNG_SIZE); std::copy(src->mm.begin(), src->mm.end(), dst.mm); dst.rngblocks = src->rngblocks; } void MCell4Converter::convert_checkpointed_molecules() { // single partition for now Partition& p = world->get_partition(PARTITION_ID_INITIAL); uint_set used_mol_ids; // add each mol for (const shared_ptr& m: model->checkpointed_molecules) { MCell::Molecule res_m; // base data if (used_mol_ids.count(m->id) != 0) { throw RuntimeError("Checkpointed molecule with ID " + to_string(m->id) + " was already added."); } res_m.id = m->id; if (m->species->species_id == SPECIES_ID_INVALID) { throw RuntimeError("Species " + m->species->name + " for checkpointed molecule is not present in the model."); } res_m.species_id = m->species->species_id; // TODO: check that the times are not in the past res_m.diffusion_time = m->diffusion_time / world->config.time_unit; res_m.birthday = m->birthday / world->config.time_unit; // TODO: check that the flags make sense res_m.flags = m->flags; if (is_set(m->unimol_rxn_time)) { res_m.unimol_rxn_time = m->unimol_rxn_time / world->config.time_unit; } else { res_m.unimol_rxn_time = TIME_INVALID; } if (m->type == MoleculeType::VOLUME) { res_m.reset_vol_data(); const shared_ptr& vm = dynamic_pointer_cast(m); res_m.v.pos = vm->pos * Vec3(world->config.rcp_length_unit); p.add_volume_molecule(res_m, 0); } else { res_m.reset_surf_data(); assert(m->type == MoleculeType::SURFACE); const shared_ptr& sm = dynamic_pointer_cast(m); res_m.s.pos = sm->pos * Vec2(world->config.rcp_length_unit); res_m.s.orientation = convert_api_orientation(sm->orientation, false); res_m.s.wall_index = sm->geometry_object->get_partition_wall_index(sm->wall_index); res_m.s.grid_tile_index = sm->grid_tile_index; // we must initialize grid and register the molecule there MCell::Wall& w = p.get_wall(res_m.s.wall_index); if (!w.has_initialized_grid()) { w.initialize_grid(p); } w.grid.set_molecule_tile(res_m.s.grid_tile_index, res_m.id); p.add_surface_molecule(res_m, 0); } } } } // namespace API } // namespace MCell