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) << 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           3 :     const auto& distribution = membraneType.discTypeDistribution;
     121             : 
     122           3 :     if (membraneType.discCount < 0)
     123           0 :         throw ExceptionWithLocation("Disc count for membrane type \"" + membraneTypeName + "\" is negative (" +
     124           0 :                                     std::to_string(membraneType.discCount) + ")");
     125             : 
     126           3 :     if (membraneType.discCount == 0 || distribution.empty())
     127           0 :         return;
     128             : 
     129           3 :     auto gridPoints = calculateCompartmentGridPoints(compartment, maxRadius, membraneType.discCount);
     130           3 :     const auto discCount = std::min(static_cast<std::size_t>(membraneType.discCount), gridPoints.size());
     131             : 
     132           3 :     if (!distribution.empty())
     133             :     {
     134           3 :         if (double sum = calculateValueSum(distribution) * 100; std::abs(sum - 100) > 1e-1)
     135           0 :             throw ExceptionWithLocation("Distribution for membrane type \"" + membraneTypeName +
     136           0 :                                         "\" doesn't add up to 100%, it adds up to " + std::to_string(sum) + "%");
     137             :     }
     138             : 
     139           9 :     for (const auto& [discTypeName, frequency] : distribution)
     140             :     {
     141           6 :         const auto count = static_cast<int>(std::round(frequency * static_cast<double>(discCount)));
     142           6 :         const auto discTypeID = simulationContext_.discTypeRegistry.getIDFor(discTypeName);
     143           6 :         const auto& discType = simulationContext_.discTypeRegistry.getByID(discTypeID);
     144             : 
     145         306 :         for (int i = 0; i < count && !gridPoints.empty(); ++i)
     146             :         {
     147         300 :             Disc newDisc(discTypeID);
     148         300 :             newDisc.setPosition(gridPoints.back());
     149         300 :             newDisc.setVelocity(
     150         300 :                 sampleVelocityFromDistribution(simulationConfig_.mostProbableSpeed, discType.getMass()));
     151             : 
     152         300 :             compartment.addDisc(std::move(newDisc));
     153         300 :             gridPoints.pop_back();
     154             :         }
     155             :     }
     156             : 
     157           5 :     for (auto& subCompartment : compartment.getCompartments())
     158           2 :         populateCompartmentWithDistribution(*subCompartment, maxRadius);
     159           3 : }
     160             : 
     161         300 : Vector2d CellPopulator::sampleVelocityFromDistribution(double mostProbableSpeed, double m) const
     162             : {
     163         300 :     static thread_local std::mt19937 gen(std::random_device{}());
     164         300 :     std::normal_distribution<double> normalDistribution(0, mostProbableSpeed / std::sqrt(m));
     165         300 :     return Vector2d{normalDistribution(gen), normalDistribution(gen)};
     166             : }
     167             : 
     168          23 : Compartment& CellPopulator::findDeepestContainingCompartment(const Disc& disc)
     169             : {
     170          23 :     const auto& M = disc.getPosition();
     171          23 :     const auto& R = simulationContext_.discTypeRegistry.getByID(disc.getTypeID()).getRadius();
     172             : 
     173          39 :     auto isFullyContainedIn = [&](const Compartment& compartment)
     174             :     {
     175          39 :         const auto& membrane = compartment.getMembrane();
     176          39 :         const auto& membraneType = simulationContext_.membraneTypeRegistry.getByID(membrane.getTypeID());
     177             : 
     178          39 :         return mathutils::circleIsFullyContainedByCircle(M, R, membrane.getPosition(), membraneType.getRadius());
     179          23 :     };
     180             : 
     181          23 :     if (!isFullyContainedIn(cell_))
     182             :     {
     183           0 :         LOG(WARNING) << "Disc at (" + std::to_string(M.x) + ", " + std::to_string(M.y) +
     184           0 :                             ") is not contained by the cell";
     185           0 :         return cell_;
     186             :     }
     187             : 
     188          23 :     Compartment* compartment = &cell_;
     189          23 :     bool continueSearch = true;
     190             : 
     191          52 :     while (continueSearch)
     192             :     {
     193          29 :         continueSearch = false;
     194          45 :         for (auto& subCompartment : compartment->getCompartments())
     195             :         {
     196          16 :             if (isFullyContainedIn(*subCompartment))
     197             :             {
     198           6 :                 compartment = subCompartment.get();
     199           6 :                 continueSearch = true;
     200             :             }
     201             :         }
     202             :     }
     203             : 
     204          23 :     return *compartment;
     205             : }
     206             : 
     207           3 : double CellPopulator::calculateValueSum(const std::unordered_map<std::string, double>& distribution) const
     208             : {
     209           3 :     return std::accumulate(distribution.begin(), distribution.end(), 0.0,
     210           9 :                            [](double partialSum, const auto& entry) { return partialSum + entry.second; });
     211             : }
     212             : 
     213             : } // namespace cell