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