Line data Source code
1 : #include "SimulationFactory.hpp" 2 : #include "Cell.hpp" 3 : #include "CellPopulator.hpp" 4 : #include "CollisionDetector.hpp" 5 : #include "CollisionHandler.hpp" 6 : #include "Disc.hpp" 7 : #include "Membrane.hpp" 8 : #include "ReactionEngine.hpp" 9 : #include "ReactionTable.hpp" 10 : #include "SimulationContext.hpp" 11 : #include "StringUtils.hpp" 12 : 13 : #include <random> 14 : 15 : namespace cell 16 : { 17 : 18 12 : SimulationFactory::SimulationFactory() = default; 19 12 : SimulationFactory::~SimulationFactory() = default; 20 : 21 13 : void SimulationFactory::buildSimulationFromConfig(const SimulationConfig& simulationConfig) 22 : { 23 : // Building might fail and we don't want anything dangling 24 13 : reset(); 25 : 26 : try 27 : { 28 13 : discTypeRegistry_ = std::make_unique<DiscTypeRegistry>(buildDiscTypeRegistry(simulationConfig)); 29 13 : membraneTypeRegistry_ = std::make_unique<MembraneTypeRegistry>(buildMembraneTypeRegistry(simulationConfig)); 30 : } 31 0 : catch (const std::exception& e) 32 : { 33 0 : throw InvalidTypesException(e.what()); 34 0 : } 35 : 36 : try 37 : { 38 : reactionTable_ = 39 13 : std::make_unique<ReactionTable>(buildReactionTable(simulationConfig, std::as_const(*discTypeRegistry_))); 40 : } 41 1 : catch (const std::exception& e) 42 : { 43 1 : throw InvalidReactionsException(e.what()); 44 1 : } 45 : 46 : try 47 : { 48 : reactionEngine_ = 49 12 : std::make_unique<ReactionEngine>(std::as_const(*discTypeRegistry_), std::as_const(*reactionTable_)); 50 24 : collisionHandler_ = std::make_unique<CollisionHandler>(std::as_const(*discTypeRegistry_), 51 24 : std::as_const(*membraneTypeRegistry_)); 52 : 53 12 : cell_ = buildCell(simulationConfig); 54 : } 55 2 : catch (const std::exception& e) 56 : { 57 2 : throw InvalidSetupException(e.what()); 58 2 : } 59 10 : } 60 : 61 45 : SimulationContext SimulationFactory::getSimulationContext() const 62 : { 63 45 : if (!discTypeRegistry_ || !membraneTypeRegistry_ || !reactionEngine_ || !collisionHandler_) 64 0 : throw ExceptionWithLocation("Can't get simulation context, dependencies haven't been fully created yet"); 65 : 66 45 : return SimulationContext{.discTypeRegistry = *discTypeRegistry_, 67 45 : .membraneTypeRegistry = *membraneTypeRegistry_, 68 45 : .reactionEngine = *reactionEngine_, 69 45 : .collisionHandler = *collisionHandler_}; 70 : } 71 : 72 10 : Cell& SimulationFactory::getCell() 73 : { 74 10 : if (!cell_) 75 0 : throw ExceptionWithLocation("Can't access cell, it hasn't been created yet"); 76 : 77 10 : return *cell_; 78 : } 79 : 80 0 : bool SimulationFactory::cellIsBuilt() const 81 : { 82 0 : return static_cast<bool>(cell_); 83 : } 84 : 85 13 : DiscTypeRegistry SimulationFactory::buildDiscTypeRegistry(const SimulationConfig& simulationConfig) 86 : { 87 13 : DiscTypeRegistry discTypeRegistry; 88 13 : std::vector<DiscType> discTypes; 89 58 : for (const auto& discType : simulationConfig.discTypes) 90 45 : discTypes.emplace_back(discType.name, Radius{discType.radius}, Mass{discType.mass}); 91 : 92 13 : discTypeRegistry.setValues(std::move(discTypes)); 93 : 94 26 : return discTypeRegistry; 95 13 : } 96 : 97 13 : MembraneTypeRegistry SimulationFactory::buildMembraneTypeRegistry(const SimulationConfig& simulationConfig) 98 : { 99 13 : MembraneTypeRegistry registry; 100 13 : std::vector<MembraneType> types; 101 : 102 26 : auto convertPermeabilityMap = [&](const std::unordered_map<std::string, MembraneType::Permeability>& configMap) 103 : { 104 26 : MembraneType::PermeabilityMap permeabilityMap; 105 33 : for (const auto& [discTypeName, permeability] : configMap) 106 7 : permeabilityMap[discTypeRegistry_->getIDFor(discTypeName)] = permeability; 107 : 108 26 : return permeabilityMap; 109 0 : }; 110 : 111 13 : const auto& cellMembraneType = simulationConfig.cellMembraneType; 112 13 : types.emplace_back(cellMembraneType.name, cellMembraneType.radius, 113 26 : convertPermeabilityMap(cellMembraneType.permeabilityMap)); 114 : 115 26 : for (const auto& type : simulationConfig.membraneTypes) 116 13 : types.emplace_back(type.name, type.radius, convertPermeabilityMap(type.permeabilityMap)); 117 : 118 13 : registry.setValues(std::move(types)); 119 : 120 26 : return registry; 121 13 : } 122 : 123 13 : ReactionTable SimulationFactory::buildReactionTable(const SimulationConfig& simulationConfig, 124 : const DiscTypeRegistry& discTypeRegistry) 125 : { 126 13 : ReactionTable reactionTable(discTypeRegistry); 127 : 128 21 : for (const auto& reaction : simulationConfig.reactions) 129 : { 130 9 : DiscTypeID educt1 = discTypeRegistry.getIDFor(reaction.educt1); 131 : std::optional<DiscTypeID> educt2 = 132 9 : reaction.educt2.empty() ? std::nullopt : std::make_optional(discTypeRegistry.getIDFor(reaction.educt2)); 133 : 134 9 : DiscTypeID product1 = discTypeRegistry.getIDFor(reaction.product1); 135 : std::optional<DiscTypeID> product2 = 136 9 : reaction.product2.empty() ? std::nullopt : std::make_optional(discTypeRegistry.getIDFor(reaction.product2)); 137 : 138 9 : Reaction newReaction(educt1, educt2, product1, product2, reaction.probability); 139 9 : if (simulationConfig.reactionsConserveArea) 140 1 : newReaction.validateAreaConservation(discTypeRegistry); 141 : 142 8 : reactionTable.addReaction(newReaction); 143 : } 144 : 145 12 : return reactionTable; 146 1 : } 147 : 148 12 : std::unique_ptr<Cell> SimulationFactory::buildCell(const SimulationConfig& simulationConfig) 149 : { 150 12 : const auto& configMembranes = simulationConfig.membranes; 151 12 : if (std::find_if(configMembranes.begin(), configMembranes.end(), [&](const auto& membrane) 152 38 : { return membrane.membraneTypeName == config::cellMembraneTypeName; }) != configMembranes.end()) 153 0 : throw ExceptionWithLocation("There can't be more than 1 cell membrane in the simulation"); 154 : 155 12 : for (const auto& [discTypeID, permeability] : simulationConfig.cellMembraneType.permeabilityMap) 156 : { 157 0 : if (permeability != MembraneType::Permeability::None) 158 0 : throw ExceptionWithLocation("Currently the outer cell membrane does not support permeability"); 159 : } 160 : 161 12 : throwIfDiscsCanBeLargerThanMembranes(simulationConfig); 162 : 163 11 : std::vector<Membrane> membranes = getMembranesFromConfig(simulationConfig); 164 : 165 11 : Membrane cellMembrane(membraneTypeRegistry_->getIDFor(config::cellMembraneTypeName)); 166 11 : cellMembrane.setPosition({0, 0}); 167 : 168 11 : std::unique_ptr<Cell> cell(std::make_unique<Cell>(std::move(cellMembrane), getSimulationContext())); 169 12 : createCompartments(*cell, std::move(membranes)); 170 : 171 10 : CellPopulator cellPopulator(*cell, simulationConfig, *discTypeRegistry_, *membraneTypeRegistry_); 172 10 : cellPopulator.populateCell(); 173 : 174 20 : return cell; 175 12 : } 176 : 177 11 : std::vector<Membrane> SimulationFactory::getMembranesFromConfig(const SimulationConfig& simulationConfig_) 178 : { 179 11 : std::vector<Membrane> membranes; 180 : 181 22 : for (const auto& membrane : simulationConfig_.membranes) 182 : { 183 11 : Membrane newMembrane(membraneTypeRegistry_->getIDFor(membrane.membraneTypeName)); 184 11 : newMembrane.setPosition({membrane.x, membrane.y}); 185 : 186 11 : membranes.push_back(std::move(newMembrane)); 187 : } 188 : 189 11 : return membranes; 190 0 : } 191 : 192 13 : void SimulationFactory::reset() 193 : { 194 13 : discTypeRegistry_.reset(); 195 13 : membraneTypeRegistry_.reset(); 196 13 : reactionTable_.reset(); 197 13 : reactionEngine_.reset(); 198 13 : collisionHandler_.reset(); 199 13 : cell_.reset(); 200 13 : } 201 : 202 11 : void SimulationFactory::createCompartments(Cell& cell, std::vector<Membrane> membranes) 203 : { 204 : // We'll sort ascending by size and convert membranes to compartments from large to small 205 : // This way, compartments are sorted descending by size, and finding a parent compartment is simple: 206 : // Iterate from small to large and find the first compartment that contains the current fully 207 : 208 11 : std::sort(membranes.begin(), membranes.end(), 209 12 : [&](const Membrane& lhs, const Membrane& rhs) 210 : { 211 12 : return membraneTypeRegistry_->getByID(lhs.getTypeID()).getRadius() < 212 12 : membraneTypeRegistry_->getByID(rhs.getTypeID()).getRadius(); 213 : }); 214 : 215 11 : std::vector<Compartment*> compartments({&cell}); 216 22 : while (!membranes.empty()) 217 : { 218 11 : auto membrane = std::move(membranes.back()); 219 11 : membranes.pop_back(); 220 : 221 11 : const auto& M = membrane.getPosition(); 222 11 : const auto& R = membraneTypeRegistry_->getByID(membrane.getTypeID()).getRadius(); 223 : 224 11 : bool parentFound = false; 225 12 : for (auto iter = compartments.rbegin(); !parentFound && iter != compartments.rend(); ++iter) 226 : { 227 12 : const auto& Mo = (*iter)->getMembrane().getPosition(); 228 12 : const auto& Ro = membraneTypeRegistry_->getByID((*iter)->getMembrane().getTypeID()).getRadius(); 229 : 230 12 : if (!mathutils::circleIsFullyContainedByCircle(M, R, Mo, Ro)) 231 1 : continue; 232 : 233 11 : auto compartment = (*iter)->createSubCompartment(std::move(membrane)); 234 11 : compartments.push_back(compartment); 235 11 : parentFound = true; 236 11 : break; 237 : } 238 : 239 11 : if (!parentFound) 240 0 : throw ExceptionWithLocation("Couldn't find a compartment that fully contains the membrane at (" + 241 0 : std::to_string(M.x) + "," + std::to_string(M.y) + ")"); 242 : } 243 : 244 : // We do it here to include the cell itself in the check 245 11 : throwIfCompartmentsIntersect(compartments); 246 11 : } 247 : 248 11 : void SimulationFactory::throwIfCompartmentsIntersect(const std::vector<Compartment*>& compartments) const 249 : { 250 : // This is done once, so no worries about O(n^2) 251 30 : for (std::size_t i = 0; i < compartments.size(); ++i) 252 : { 253 20 : const auto R1 = membraneTypeRegistry_->getByID(compartments[i]->getMembrane().getTypeID()).getRadius(); 254 20 : const auto M1 = compartments[i]->getMembrane().getPosition(); 255 : 256 36 : for (std::size_t j = i + 1; j < compartments.size(); ++j) 257 : { 258 17 : const auto R2 = membraneTypeRegistry_->getByID(compartments[j]->getMembrane().getTypeID()).getRadius(); 259 17 : const auto M2 = compartments[j]->getMembrane().getPosition(); 260 : 261 17 : if (mathutils::circlesIntersect(M1, R1, M2, R2)) 262 3 : throw ExceptionWithLocation("Overlapping membranes at " + stringutils::toString(M1) + " and " + 263 4 : stringutils::toString(M2) + " not allowed"); 264 : } 265 : } 266 10 : } 267 : 268 12 : void SimulationFactory::throwIfDiscsCanBeLargerThanMembranes(const SimulationConfig& config) const 269 : { 270 12 : if (config.discTypes.empty()) 271 0 : return; 272 : 273 12 : const auto maxDiscRadius = std::max_element(config.discTypes.begin(), config.discTypes.end(), 274 29 : [](const auto& t1, const auto& t2) { return t1.radius < t2.radius; }) 275 12 : ->radius; 276 : 277 12 : auto minMembraneRadius = config.cellMembraneType.radius; 278 25 : for (const auto& membraneType : config.membraneTypes) 279 13 : minMembraneRadius = std::min(minMembraneRadius, membraneType.radius); 280 : 281 12 : if (maxDiscRadius >= minMembraneRadius) 282 1 : throw ExceptionWithLocation("At least one disc type has a radius larger than at least one membrane type."); 283 : } 284 : 285 : } // namespace cell