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