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