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           9 : CellPopulator::CellPopulator(Cell& cell, SimulationConfig simulationConfig, SimulationContext simulationContext)
      17           9 :     : cell_(cell)
      18           9 :     , simulationConfig_(std::move(simulationConfig))
      19           9 :     , simulationContext_(std::move(simulationContext))
      20             : {
      21           9 : }
      22             : 
      23           9 : void CellPopulator::populateCell()
      24             : {
      25           9 :     if (simulationConfig_.useDistribution)
      26           1 :         populateWithDistributions();
      27             :     else
      28           8 :         populateDirectly();
      29           9 : }
      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           8 : void CellPopulator::populateDirectly()
      46             : {
      47          31 :     for (const auto& disc : simulationConfig_.discs)
      48             :     {
      49          23 :         Disc newDisc(simulationContext_.discTypeRegistry.getIDFor(disc.discTypeName));
      50          23 :         newDisc.setPosition({disc.x, disc.y});
      51          23 :         newDisc.setVelocity({disc.vx, disc.vy});
      52             : 
      53          23 :         auto& compartment = findDeepestContainingCompartment(newDisc);
      54          23 :         compartment.addDisc(std::move(newDisc));
      55             :     }
      56           8 : }
      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<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 - 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) << discCount << " discs should be created for membrane of type \"" << membraneType.getName()
     105           0 :                      << "\", but the grid can only fit " << gridPoints.size() << ". " << discCount - gridPoints.size()
     106           0 :                      << " discs will not be created.";
     107             :     }
     108             : 
     109           6 :     return gridPoints;
     110           0 : }
     111             : 
     112           3 : void CellPopulator::populateCompartmentWithDistribution(Compartment& compartment, double maxRadius)
     113             : {
     114           5 :     for (auto& subCompartment : compartment.getCompartments())
     115           2 :         populateCompartmentWithDistribution(*subCompartment, maxRadius);
     116             : 
     117             :     const auto& membraneTypeName =
     118           3 :         simulationContext_.membraneTypeRegistry.getByID(compartment.getMembrane().getTypeID()).getName();
     119             : 
     120           3 :     const auto& membraneType = findMembraneTypeByName(simulationConfig_, membraneTypeName);
     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 :     if (membraneType.discCount == 0 || distribution.empty())
     128           0 :         return;
     129             : 
     130           3 :     auto gridPoints = calculateCompartmentGridPoints(compartment, maxRadius, membraneType.discCount);
     131           3 :     const auto discCount = std::min(static_cast<std::size_t>(membraneType.discCount), gridPoints.size());
     132             : 
     133           3 :     if (!distribution.empty())
     134             :     {
     135           3 :         if (double sum = calculateValueSum(distribution) * 100; std::abs(sum - 100) > 1e-1)
     136           0 :             throw ExceptionWithLocation("Distribution for membrane type \"" + membraneTypeName +
     137           0 :                                         "\" doesn't add up to 100%, it adds up to " + std::to_string(sum) + "%");
     138             :     }
     139             : 
     140           9 :     for (const auto& [discTypeName, frequency] : distribution)
     141             :     {
     142           6 :         const auto count = static_cast<int>(std::round(frequency * static_cast<double>(discCount)));
     143           6 :         const auto discTypeID = simulationContext_.discTypeRegistry.getIDFor(discTypeName);
     144           6 :         const auto& discType = simulationContext_.discTypeRegistry.getByID(discTypeID);
     145             : 
     146         306 :         for (int i = 0; i < count && !gridPoints.empty(); ++i)
     147             :         {
     148         300 :             Disc newDisc(discTypeID);
     149         300 :             newDisc.setPosition(gridPoints.back());
     150         300 :             newDisc.setVelocity(
     151         300 :                 sampleVelocityFromDistribution(simulationConfig_.mostProbableSpeed, discType.getMass()));
     152             : 
     153         300 :             compartment.addDisc(std::move(newDisc));
     154         300 :             gridPoints.pop_back();
     155             :         }
     156             :     }
     157           3 : }
     158             : 
     159         300 : Vector2d CellPopulator::sampleVelocityFromDistribution(double mostProbableSpeed, double m) const
     160             : {
     161         300 :     static thread_local std::mt19937 gen(std::random_device{}());
     162         300 :     std::normal_distribution<double> normalDistribution(0, mostProbableSpeed / std::sqrt(m));
     163         300 :     return Vector2d{normalDistribution(gen), normalDistribution(gen)};
     164             : }
     165             : 
     166          23 : Compartment& CellPopulator::findDeepestContainingCompartment(const Disc& disc)
     167             : {
     168          23 :     const auto& M = disc.getPosition();
     169          23 :     const auto& R = simulationContext_.discTypeRegistry.getByID(disc.getTypeID()).getRadius();
     170             : 
     171          39 :     auto isFullyContainedIn = [&](const Compartment& compartment)
     172             :     {
     173          39 :         const auto& membrane = compartment.getMembrane();
     174          39 :         const auto& membraneType = simulationContext_.membraneTypeRegistry.getByID(membrane.getTypeID());
     175             : 
     176          39 :         return mathutils::circleIsFullyContainedByCircle(M, R, membrane.getPosition(), membraneType.getRadius());
     177          23 :     };
     178             : 
     179          23 :     if (!isFullyContainedIn(cell_))
     180             :     {
     181           0 :         LOG(WARNING) << "Disc at (" + std::to_string(M.x) + ", " + std::to_string(M.y) +
     182           0 :                             ") is not contained by the cell";
     183           0 :         return cell_;
     184             :     }
     185             : 
     186          23 :     Compartment* compartment = &cell_;
     187          23 :     bool continueSearch = true;
     188             : 
     189          52 :     while (continueSearch)
     190             :     {
     191          29 :         continueSearch = false;
     192          45 :         for (auto& subCompartment : compartment->getCompartments())
     193             :         {
     194          16 :             if (isFullyContainedIn(*subCompartment))
     195             :             {
     196           6 :                 compartment = subCompartment.get();
     197           6 :                 continueSearch = true;
     198             :             }
     199             :         }
     200             :     }
     201             : 
     202          23 :     return *compartment;
     203             : }
     204             : 
     205           3 : double CellPopulator::calculateValueSum(const std::unordered_map<std::string, double>& distribution) const
     206             : {
     207           3 :     return std::accumulate(distribution.begin(), distribution.end(), 0.0,
     208           9 :                            [](double partialSum, const auto& entry) { return partialSum + entry.second; });
     209             : }
     210             : 
     211             : } // namespace cell