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