Line data    Source code

       1             : #include "CellPopulator.hpp"
       2             : #include "Cell.hpp"
       3             : #include "Disc.hpp"
       4             : #include "MathUtils.hpp"
       5             : 
       6             : #include <glog/logging.h>
       7             : 
       8             : #include <algorithm>
       9             : #include <numbers>
      10             : #include <random>
      11             : #include <vector>
      12             : 
      13             : namespace cell
      14             : {
      15             : 
      16           8 : CellPopulator::CellPopulator(Cell& cell, SimulationConfig simulationConfig, SimulationContext simulationContext)
      17           8 :     : cell_(cell)
      18           8 :     , simulationConfig_(std::move(simulationConfig))
      19           8 :     , simulationContext_(std::move(simulationContext))
      20             : {
      21           8 : }
      22             : 
      23           8 : void CellPopulator::populateCell()
      24             : {
      25           8 :     if (simulationConfig_.useDistribution)
      26           1 :         populateWithDistributions();
      27             :     else
      28           7 :         populateDirectly();
      29           8 : }
      30             : 
      31           1 : void CellPopulator::populateWithDistributions()
      32             : {
      33           1 :     const auto& discTypes = simulationConfig_.discTypes;
      34           1 :     if (discTypes.empty())
      35           0 :         return;
      36             : 
      37           1 :     double maxRadius = std::max_element(discTypes.begin(), discTypes.end(),
      38           1 :                                         [](const config::DiscType& lhs, const config::DiscType& rhs)
      39           1 :                                         { return lhs.radius < rhs.radius; })
      40           1 :                            ->radius;
      41             : 
      42           1 :     populateCompartmentWithDistribution(cell_, maxRadius);
      43             : }
      44             : 
      45           7 : void CellPopulator::populateDirectly()
      46             : {
      47          26 :     for (const auto& disc : simulationConfig_.discs)
      48             :     {
      49          19 :         Disc newDisc(simulationContext_.discTypeRegistry.getIDFor(disc.discTypeName));
      50          19 :         newDisc.setPosition({disc.x, disc.y});
      51          19 :         newDisc.setVelocity({disc.vx, disc.vy});
      52             : 
      53          19 :         auto& compartment = findDeepestContainingCompartment(newDisc);
      54          19 :         compartment.addDisc(std::move(newDisc));
      55             :     }
      56           7 : }
      57             : 
      58           0 : double CellPopulator::calculateDistributionSum(const std::map<std::string, double>& distribution) const
      59             : {
      60           0 :     return std::accumulate(distribution.begin(), distribution.end(), 0.0,
      61           0 :                            [](double currentSum, auto& entryPair) { return currentSum + entryPair.second; });
      62             : }
      63             : 
      64           3 : std::vector<sf::Vector2d> CellPopulator::calculateCompartmentGridPoints(Compartment& compartment, double maxRadius,
      65             :                                                                         int discCount) const
      66             : {
      67           3 :     const auto& membraneType = simulationContext_.membraneTypeRegistry.getByID(compartment.getMembrane().getTypeID());
      68           3 :     const auto& membraneCenter = compartment.getMembrane().getPosition();
      69           3 :     const auto& membraneRadius = membraneType.getRadius();
      70             : 
      71           3 :     auto gridPoints = mathutils::calculateGrid(2 * membraneRadius, 2 * membraneRadius, 2 * maxRadius);
      72             : 
      73           3 :     const auto& topLeft = membraneCenter - sf::Vector2d{membraneRadius, membraneRadius};
      74      462337 :     for (size_t i = 0; i < gridPoints.size();)
      75             :     {
      76      462334 :         gridPoints[i] += topLeft;
      77      462334 :         bool valid = true;
      78             : 
      79      462334 :         if (!mathutils::circleIsFullyContainedByCircle(gridPoints[i], maxRadius, membraneCenter, membraneRadius))
      80       99320 :             valid = false;
      81             : 
      82      462334 :         const auto& subCompartments = compartment.getCompartments();
      83      824515 :         for (auto iter = subCompartments.begin(); valid && iter != subCompartments.end(); ++iter)
      84             :         {
      85      362181 :             const auto& membrane = (*iter)->getMembrane();
      86      362181 :             const auto& M = membrane.getPosition();
      87      362181 :             const auto& R = simulationContext_.membraneTypeRegistry.getByID(membrane.getTypeID()).getRadius();
      88             : 
      89      362181 :             if (mathutils::circlesOverlap(gridPoints[i], maxRadius, M, R))
      90       14985 :                 valid = false;
      91             :         }
      92             : 
      93      462334 :         if (!valid)
      94             :         {
      95      114305 :             std::swap(gridPoints[i], gridPoints.back());
      96      114305 :             gridPoints.pop_back();
      97             :         }
      98             :         else
      99      348029 :             ++i;
     100             :     }
     101             : 
     102           3 :     if (static_cast<int>(gridPoints.size()) < discCount)
     103             :     {
     104           0 :         LOG(WARNING) << std::to_string(discCount)
     105           0 :                      << " discs should be created for membrane of type \"" + membraneType.getName() +
     106             :                             "\", but the grid can only fit "
     107           0 :                      << std::to_string(gridPoints.size()) << ". " << std::to_string(discCount - gridPoints.size())
     108           0 :                      << " discs will not be created.";
     109             :     }
     110             : 
     111           6 :     return gridPoints;
     112           0 : }
     113             : 
     114           3 : void CellPopulator::populateCompartmentWithDistribution(Compartment& compartment, double maxRadius)
     115             : {
     116             :     const auto& membraneTypeName =
     117           3 :         simulationContext_.membraneTypeRegistry.getByID(compartment.getMembrane().getTypeID()).getName();
     118             : 
     119           3 :     const auto& membraneType = findMembraneTypeByName(simulationConfig_, membraneTypeName);
     120             : 
     121           3 :     const auto& distribution = membraneType.discTypeDistribution;
     122             : 
     123           3 :     if (membraneType.discCount < 0)
     124           0 :         throw ExceptionWithLocation("Disc count for membrane type \"" + membraneTypeName + "\" is negative (" +
     125           0 :                                     std::to_string(membraneType.discCount) + ")");
     126             : 
     127           3 :     auto gridPoints = calculateCompartmentGridPoints(compartment, maxRadius, membraneType.discCount);
     128           3 :     const auto discCount = std::min(static_cast<std::size_t>(membraneType.discCount), gridPoints.size());
     129             : 
     130           3 :     if (!distribution.empty())
     131             :     {
     132           3 :         if (double sum = calculateValueSum(distribution) * 100; std::abs(sum - 100) > 1e-1)
     133           0 :             throw ExceptionWithLocation("Distribution for membrane type \"" + membraneTypeName +
     134           0 :                                         "\" doesn't add up to 100%, it adds up to " + std::to_string(sum) + "%");
     135             :     }
     136             : 
     137           9 :     for (const auto& [discTypeName, frequency] : distribution)
     138             :     {
     139           6 :         const auto count = static_cast<int>(std::round(frequency * static_cast<double>(discCount)));
     140           6 :         const auto discTypeID = simulationContext_.discTypeRegistry.getIDFor(discTypeName);
     141             : 
     142         306 :         for (int i = 0; i < count && !gridPoints.empty(); ++i)
     143             :         {
     144         300 :             Disc newDisc(discTypeID);
     145         300 :             newDisc.setPosition(gridPoints.back());
     146         300 :             newDisc.setVelocity(sampleVelocityFromDistribution());
     147             : 
     148         300 :             compartment.addDisc(std::move(newDisc));
     149         300 :             gridPoints.pop_back();
     150             :         }
     151             :     }
     152             : 
     153           5 :     for (auto& subCompartment : compartment.getCompartments())
     154           2 :         populateCompartmentWithDistribution(*subCompartment, maxRadius);
     155           3 : }
     156             : 
     157         300 : sf::Vector2d CellPopulator::sampleVelocityFromDistribution() const
     158             : {
     159         300 :     static thread_local std::random_device rd;
     160         300 :     static thread_local std::mt19937 gen(rd());
     161             : 
     162         300 :     const auto sigma = simulationConfig_.maxVelocity;
     163             : 
     164             :     // 2D maxwell distribution is a rayleigh distribution
     165             :     // weibull distribution with k=2 and lambda = sigma*sqrt(2) is rayleigh distribution
     166             : 
     167         300 :     std::uniform_real_distribution<double> angleDistribution(0, 2 * std::numbers::pi);
     168         300 :     std::weibull_distribution<double> speedDistribution(2, std::numbers::sqrt2 * sigma);
     169             : 
     170         300 :     const auto angle = angleDistribution(gen);
     171         300 :     const auto speed = speedDistribution(gen);
     172             : 
     173         300 :     return sf::Vector2d{speed * std::cos(angle), speed * std::sin(angle)};
     174             : }
     175             : 
     176          19 : Compartment& CellPopulator::findDeepestContainingCompartment(const Disc& disc)
     177             : {
     178          19 :     const auto& M = disc.getPosition();
     179          19 :     const auto& R = simulationContext_.discTypeRegistry.getByID(disc.getTypeID()).getRadius();
     180             : 
     181          31 :     auto isFullyContainedIn = [&](const Compartment& compartment)
     182             :     {
     183          31 :         const auto& membrane = compartment.getMembrane();
     184          31 :         const auto& membraneType = simulationContext_.membraneTypeRegistry.getByID(membrane.getTypeID());
     185             : 
     186          31 :         return mathutils::circleIsFullyContainedByCircle(M, R, membrane.getPosition(), membraneType.getRadius());
     187          19 :     };
     188             : 
     189          19 :     if (!isFullyContainedIn(cell_))
     190             :     {
     191           0 :         LOG(WARNING) << "Disc at (" + std::to_string(M.x) + ", " + std::to_string(M.y) +
     192           0 :                             ") is not contained by the cell";
     193           0 :         return cell_;
     194             :     }
     195             : 
     196          19 :     Compartment* compartment = &cell_;
     197          19 :     bool continueSearch = true;
     198             : 
     199          42 :     while (continueSearch)
     200             :     {
     201          23 :         continueSearch = false;
     202          35 :         for (auto& subCompartment : compartment->getCompartments())
     203             :         {
     204          12 :             if (isFullyContainedIn(*subCompartment))
     205             :             {
     206           4 :                 compartment = subCompartment.get();
     207           4 :                 continueSearch = true;
     208             :             }
     209             :         }
     210             :     }
     211             : 
     212          19 :     return *compartment;
     213             : }
     214             : 
     215           3 : double CellPopulator::calculateValueSum(const std::unordered_map<std::string, double>& distribution) const
     216             : {
     217           3 :     return std::accumulate(distribution.begin(), distribution.end(), 0.0,
     218           9 :                            [](double partialSum, const auto& entry) { return partialSum + entry.second; });
     219             : }
     220             : 
     221             : } // namespace cell