Line data    Source code

       1             : #include "CellPopulator.hpp"
       2             : #include "Cell.hpp"
       3             : #include "Disc.hpp"
       4             : #include "MathUtils.hpp"
       5             : 
       6             : #include <algorithm>
       7             : #include <iostream>
       8             : #include <numbers>
       9             : #include <random>
      10             : #include <vector>
      11             : 
      12             : namespace cell
      13             : {
      14             : 
      15           9 : CellPopulator::CellPopulator(Cell& cell, SimulationConfig simulationConfig, const DiscTypeRegistry& discTypeRegistry,
      16           9 :                              const MembraneTypeRegistry& membraneTypeRegistry)
      17           9 :     : cell_(cell)
      18           9 :     , simulationConfig_(std::move(simulationConfig))
      19           9 :     , discTypeRegistry_(discTypeRegistry)
      20           9 :     , membraneTypeRegistry_(membraneTypeRegistry)
      21             : {
      22           9 : }
      23             : 
      24           9 : void CellPopulator::populateCell()
      25             : {
      26           9 :     if (simulationConfig_.useDistribution)
      27           1 :         populateWithDistributions();
      28             :     else
      29           8 :         populateDirectly();
      30           9 : }
      31             : 
      32           1 : void CellPopulator::populateWithDistributions()
      33             : {
      34           1 :     const auto& discTypes = simulationConfig_.discTypes;
      35           1 :     if (discTypes.empty())
      36           0 :         return;
      37             : 
      38           1 :     double maxRadius = std::max_element(discTypes.begin(), discTypes.end(),
      39           1 :                                         [](const config::DiscType& lhs, const config::DiscType& rhs)
      40           1 :                                         { return lhs.radius < rhs.radius; })
      41           1 :                            ->radius;
      42             : 
      43           1 :     populateCompartmentWithDistribution(cell_, maxRadius);
      44             : }
      45             : 
      46           8 : void CellPopulator::populateDirectly()
      47             : {
      48          31 :     for (const auto& disc : simulationConfig_.discs)
      49             :     {
      50          23 :         Disc newDisc(discTypeRegistry_.getIDFor(disc.discTypeName));
      51          23 :         newDisc.setPosition({disc.x, disc.y});
      52          23 :         newDisc.setVelocity({disc.vx, disc.vy});
      53             : 
      54          23 :         auto& compartment = findDeepestContainingCompartment(newDisc);
      55          23 :         compartment.addDisc(std::move(newDisc));
      56             :     }
      57           8 : }
      58             : 
      59           0 : double CellPopulator::calculateDistributionSum(const std::map<std::string, double>& distribution) const
      60             : {
      61           0 :     return std::accumulate(distribution.begin(), distribution.end(), 0.0,
      62           0 :                            [](double currentSum, auto& entryPair) { return currentSum + entryPair.second; });
      63             : }
      64             : 
      65           3 : std::vector<Vector2d> CellPopulator::calculateCompartmentGridPoints(Compartment& compartment, double maxRadius,
      66             :                                                                     int discCount) const
      67             : {
      68           3 :     const auto& membraneType = membraneTypeRegistry_.getByID(compartment.getMembrane().getTypeID());
      69           3 :     const auto& membraneCenter = compartment.getMembrane().getPosition();
      70           3 :     const auto& membraneRadius = membraneType.getRadius();
      71             : 
      72           3 :     auto gridPoints = mathutils::calculateGrid(2 * membraneRadius, 2 * membraneRadius, 2 * maxRadius);
      73             : 
      74           3 :     const auto& topLeft = membraneCenter - Vector2d{membraneRadius, membraneRadius};
      75      462337 :     for (size_t i = 0; i < gridPoints.size();)
      76             :     {
      77      462334 :         gridPoints[i] += topLeft;
      78      462334 :         bool valid = true;
      79             : 
      80      462334 :         if (!mathutils::circleIsFullyContainedByCircle(gridPoints[i], maxRadius, membraneCenter, membraneRadius))
      81       99320 :             valid = false;
      82             : 
      83      462334 :         const auto& subCompartments = compartment.getCompartments();
      84      824515 :         for (auto iter = subCompartments.begin(); valid && iter != subCompartments.end(); ++iter)
      85             :         {
      86      362181 :             const auto& membrane = (*iter)->getMembrane();
      87      362181 :             const auto& M = membrane.getPosition();
      88      362181 :             const auto& R = membraneTypeRegistry_.getByID(membrane.getTypeID()).getRadius();
      89             : 
      90      362181 :             if (mathutils::circlesOverlap(gridPoints[i], maxRadius, M, R))
      91       14985 :                 valid = false;
      92             :         }
      93             : 
      94      462334 :         if (!valid)
      95             :         {
      96      114305 :             std::swap(gridPoints[i], gridPoints.back());
      97      114305 :             gridPoints.pop_back();
      98             :         }
      99             :         else
     100      348029 :             ++i;
     101             :     }
     102             : 
     103           3 :     if (static_cast<int>(gridPoints.size()) < discCount)
     104             :     {
     105           0 :         std::cout << "Grid for \"" << membraneType.getName() << "\" can only fit " << gridPoints.size() << "/"
     106           0 :                   << discCount << " discs\n";
     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           3 :     const auto& membraneTypeName = 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 = discTypeRegistry_.getIDFor(discTypeName);
     143           6 :         const auto& discType = 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           3 : }
     157             : 
     158         300 : Vector2d CellPopulator::sampleVelocityFromDistribution(double mostProbableSpeed, double m) const
     159             : {
     160         300 :     static thread_local std::mt19937 gen(std::random_device{}());
     161         300 :     std::normal_distribution<double> normalDistribution(0, mostProbableSpeed / std::sqrt(m));
     162         300 :     return Vector2d{normalDistribution(gen), normalDistribution(gen)};
     163             : }
     164             : 
     165          23 : Compartment& CellPopulator::findDeepestContainingCompartment(const Disc& disc)
     166             : {
     167          23 :     const auto& M = disc.getPosition();
     168          23 :     const auto& R = discTypeRegistry_.getByID(disc.getTypeID()).getRadius();
     169             : 
     170          39 :     auto isFullyContainedIn = [&](const Compartment& compartment)
     171             :     {
     172          39 :         const auto& membrane = compartment.getMembrane();
     173          39 :         const auto& membraneType = membraneTypeRegistry_.getByID(membrane.getTypeID());
     174             : 
     175          39 :         return mathutils::circleIsFullyContainedByCircle(M, R, membrane.getPosition(), membraneType.getRadius());
     176          23 :     };
     177             : 
     178          23 :     if (!isFullyContainedIn(cell_))
     179             :     {
     180           0 :         std::cout << "Disc at (" + std::to_string(M.x) + ", " + std::to_string(M.y) + ") is not contained by the cell";
     181           0 :         return cell_;
     182             :     }
     183             : 
     184          23 :     Compartment* compartment = &cell_;
     185          23 :     bool continueSearch = true;
     186             : 
     187          52 :     while (continueSearch)
     188             :     {
     189          29 :         continueSearch = false;
     190          45 :         for (auto& subCompartment : compartment->getCompartments())
     191             :         {
     192          16 :             if (isFullyContainedIn(*subCompartment))
     193             :             {
     194           6 :                 compartment = subCompartment.get();
     195           6 :                 continueSearch = true;
     196             :             }
     197             :         }
     198             :     }
     199             : 
     200          23 :     return *compartment;
     201             : }
     202             : 
     203           3 : double CellPopulator::calculateValueSum(const std::unordered_map<std::string, double>& distribution) const
     204             : {
     205           3 :     return std::accumulate(distribution.begin(), distribution.end(), 0.0,
     206           9 :                            [](double partialSum, const auto& entry) { return partialSum + entry.second; });
     207             : }
     208             : 
     209             : } // namespace cell