Area conservation and missing histogram plots

scripts/unix/clang-format.sh

@@ -5,8 +5,8 @@ cd ../..
 
 if [[ "$1" == "-f" ]]; then
     echo "Formatting files..."
-    find src/ include/ test/ -name '*.h' -o -name '*.hpp' -o -name '*.inl' -o -name '*.cpp' | xargs clang-format -i
+    find src/ test/ -name '*.h' -o -name '*.hpp' -o -name '*.inl' -o -name '*.cpp' | xargs clang-format -i
 else
     echo "Checking with clang-format..."
-    find src/ include/ test/ -name '*.h' -o -name '*.hpp' -o -name '*.inl' -o -name '*.cpp' | xargs clang-format --dry-run --Werror
+    find src/ test/ -name '*.h' -o -name '*.hpp' -o -name '*.inl' -o -name '*.cpp' | xargs clang-format --dry-run --Werror
 fi

src/apps/cell-cli/Main.cpp

@@ -14,7 +14,7 @@ using namespace std::chrono_literals;
 
 int main(int argc, char** argv)
 {
-    CLI::App app{"Command line interface for the cell simulation"};
+    CLI::App app{"Cell 1.1.1\nCommand line interface for the cell simulation"};
 
     fs::path configFile;
     fs::path outFile;

src/lib/cell/DataPoint.cpp

@@ -60,7 +60,7 @@ void DataPoint::initializeHistograms(const std::vector<DiscTypeID>& discTypeIDs,
                                            bh::axis::regular<>(20, -3 * vSigma, 3 * vSigma, "v_y"));
 
     data_.vHistogram = bh::make_histogram(bh::axis::category<DiscTypeID>(discTypeIDs, "Disc type"),
-                                          bh::axis::regular<>(20, -3 * vSigma, 3 * vSigma, "v"));
+                                          bh::axis::regular<>(20, 0, 4 * vSigma, "v"));
 }
 
 void DataPoint::addSimulationData(Cell& cell, const ch::duration<double>& elapsedTime,

src/lib/cell/Reaction.cpp

@@ -2,6 +2,8 @@
 #include "ExceptionWithLocation.hpp"
 #include "Hashing.hpp"
 
+#include <numbers>
+
 namespace cell
 {
 
@@ -66,19 +68,33 @@ void Reaction::setProbability(double probability)
 
 void Reaction::validate(const DiscTypeRegistry& discTypeRegistry) const
 {
-    auto getMass = [&](DiscTypeID discTypeID) { return discTypeRegistry.getByID(discTypeID).getMass(); };
+    const auto getMass = [&](DiscTypeID discTypeID) { return discTypeRegistry.getByID(discTypeID).getMass(); };
 
     const auto eductMassSum = getMass(educt1_) + (educt2_ ? getMass(*educt2_) : 0);
     const auto productMassSum = getMass(product1_) + (product2_ ? getMass(*product2_) : 0);
 
-    if (eductMassSum != productMassSum)
+    if (std::abs(eductMassSum - productMassSum) > 1e-6)
         throw ExceptionWithLocation(toString(*this, discTypeRegistry) +
                                     ": Product- and educt masses need to be identical");
 
     if (type_ == Type::Transformation && educt1_ == product1_)
         throw ExceptionWithLocation(toString(*this, discTypeRegistry) + ": Educt 1 and product 1 are identical");
 }
 
+void Reaction::validateAreaConservation(const DiscTypeRegistry& discTypeRegistry) const
+{
+    const auto getArea = [&](DiscTypeID discTypeID)
+    {
+        const auto R = discTypeRegistry.getByID(discTypeID).getRadius();
+        return std::numbers::pi * R * R;
+    };
+    const auto eductAreaSum = getArea(educt1_) + (educt2_ ? getArea(*educt2_) : 0);
+    const auto productAreaSum = getArea(product1_) + (product2_ ? getArea(*product2_) : 0);
+
+    if (std::abs(eductAreaSum - productAreaSum) > 1e-6)
+        throw ExceptionWithLocation(toString(*this, discTypeRegistry) + ": Area not conserved");
+}
+
 std::string Reaction::getTypeString() const
 {
     switch (type_)

src/lib/cell/Reaction.hpp

@@ -157,6 +157,8 @@ class Reaction
      */
     void validate(const DiscTypeRegistry& discTypeRegistry) const;
 
+    void validateAreaConservation(const DiscTypeRegistry& discTypeRegistry) const;
+
 private:
     std::string getTypeString() const;
 

src/lib/cell/SimulationConfig.hpp

@@ -91,12 +91,11 @@ struct SimulationConfig
      */
     double simulationTimeScale = 1;
     double mostProbableSpeed = 600;
-
     bool useDistribution = true;
+    bool reactionsConserveArea = false;
 
     // In case of no distribution, these are used
     std::vector<config::Disc> discs;
-
     // These never use a distribution
     std::vector<config::Membrane> membranes;
 
@@ -114,9 +113,9 @@ NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(Reaction, educt1, educt2, product1, product2,
 
 } // namespace config
 
-NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE(SimulationConfig, discTypes, membraneTypes, reactions, cellMembraneType,
-                                   simulationTimeStep, simulationTimeScale, mostProbableSpeed, useDistribution, discs,
-                                   membranes)
+NLOHMANN_DEFINE_TYPE_NON_INTRUSIVE_WITH_DEFAULT(SimulationConfig, discTypes, membraneTypes, reactions, cellMembraneType,
+                                                simulationTimeStep, simulationTimeScale, mostProbableSpeed,
+                                                useDistribution, reactionsConserveArea, discs, membranes)
 
 cell::config::MembraneType& findMembraneTypeByName(cell::SimulationConfig& simulationConfig,
                                                    std::string membraneTypeName);

src/lib/cell/SimulationConfigBuilder.cpp

@@ -80,6 +80,11 @@ void SimulationConfigBuilder::setMostProbableSpeed(double mostProbableSpeed)
     simulationConfig_.mostProbableSpeed = mostProbableSpeed;
 }
 
+void SimulationConfigBuilder::setReactionsConserveArea(bool value)
+{
+    simulationConfig_.reactionsConserveArea = value;
+}
+
 const SimulationConfig& SimulationConfigBuilder::getSimulationConfig() const
 {
     return simulationConfig_;

src/lib/cell/SimulationConfigBuilder.hpp

@@ -31,6 +31,7 @@ class SimulationConfigBuilder
     void setTimeStep(double simulationTimeStep);
     void setTimeScale(double simulationTimeScale);
     void setMostProbableSpeed(double mostProbableSpeed);
+    void setReactionsConserveArea(bool value);
 
     const SimulationConfig& getSimulationConfig() const;
 

src/lib/cell/SimulationFactory.cpp

@@ -136,6 +136,9 @@ ReactionTable SimulationFactory::buildReactionTable(const SimulationConfig& simu
             reaction.product2.empty() ? std::nullopt : std::make_optional(discTypeRegistry.getIDFor(reaction.product2));
 
         Reaction newReaction(educt1, educt2, product1, product2, reaction.probability);
+        if (simulationConfig.reactionsConserveArea)
+            newReaction.validateAreaConservation(discTypeRegistry);
+
         reactionTable.addReaction(newReaction);
     }
 

src/lib/cell/SimulationRecordSerializer.cpp

@@ -6,7 +6,7 @@
 namespace cell
 {
 
-void SimulationRecordSerializer::writeTypeCountsToCsv(const std::vector<DataPoint>& dataPoints,
+void SimulationRecordSerializer::writeTypeCountsToCsv(const std::deque<DataPoint>& dataPoints,
                                                       const DiscTypeRegistry& discTypeRegistry, const fs::path& outFile)
 {
     std::ofstream file(outFile);

src/lib/cell/SimulationRecordSerializer.hpp

@@ -4,13 +4,15 @@
 #include "DataPoint.hpp"
 #include "Types.hpp"
 
+#include <deque>
+
 namespace cell
 {
 
 class SimulationRecordSerializer
 {
 public:
-    void writeTypeCountsToCsv(const std::vector<DataPoint>& dataPoints, const DiscTypeRegistry& discTypeRegistry,
+    void writeTypeCountsToCsv(const std::deque<DataPoint>& dataPoints, const DiscTypeRegistry& discTypeRegistry,
                               const fs::path& outFile);
 };
 

src/lib/cell/SimulationRecorder.cpp

@@ -11,11 +11,6 @@ namespace cell
 SimulationRecorder::SimulationRecorder(const DiscTypeRegistry& discTypeRegistry, double vSigma)
     : discTypeRegistry_(discTypeRegistry)
 {
-    // With a simulation time step of 1ms, we get 1000 data points each second
-    // With an averaging time of 10ms, we save 100 datapoints for 1 second
-    // We'll reserve enough space for 5 minutes, 5*60*100
-    dataPoints_.reserve(30000);
-
     std::vector<DiscTypeID> discTypeIDs = discTypeRegistry.getIDs();
     currentDataPoint_.initializeHistograms(discTypeIDs, vSigma);
 }
@@ -54,7 +49,7 @@ void SimulationRecorder::storeRemainingData()
     storeDataPoint();
 }
 
-const std::vector<DataPoint>& SimulationRecorder::getDataPoints() const
+const std::deque<DataPoint>& SimulationRecorder::getDataPoints() const
 {
     return dataPoints_;
 }

src/lib/cell/SimulationRecorder.hpp

@@ -8,6 +8,8 @@
 
 #include <boost/histogram.hpp>
 
+#include <deque>
+
 namespace cell
 {
 
@@ -33,7 +35,7 @@ class SimulationRecorder
     void processInitialSimulationData(Cell& cell);
     void processSimulationData(Cell& cell, const ch::duration<double>& elapsedTime);
     void storeRemainingData();
-    const std::vector<DataPoint>& getDataPoints() const;
+    const std::deque<DataPoint>& getDataPoints() const;
     void clear();
     const DataPoint& getCurrentDataPoint() const;
     void setRecordLastFrame(bool value);
@@ -48,7 +50,7 @@ class SimulationRecorder
 private:
     ch::duration<double> storageInterval_ = ch::milliseconds{100};
     DataPoint currentDataPoint_;
-    std::vector<DataPoint> dataPoints_;
+    std::deque<DataPoint> dataPoints_;
     const DiscTypeRegistry& discTypeRegistry_;
     bool recordLastFrame_ = false;
     Frame lastFrame_;

src/lib/gui/core/MainWindow.cpp

@@ -74,7 +74,7 @@ MainWindow::MainWindow(QWidget* parent)
     connect(ui->simulationControlWidget, &SimulationControlWidget::fitIntoViewRequested, ui->simulationWidget,
             &SimulationWidget::fitSimulationIntoView);
 
-    connect(simulation_.get(), &Simulation::initialFrame, ui->simulationWidget,
+    connect(simulation_.get(), &Simulation::drawFrameImmediately, ui->simulationWidget,
             &SimulationWidget::renderFrameImmediately);
     connect(simulation_.get(), &Simulation::frame, ui->simulationWidget, &SimulationWidget::queueFrameForRendering);
     connect(simulation_.get(), &Simulation::performanceData, ui->simulationInfoWidget,
@@ -192,7 +192,7 @@ void MainWindow::showAboutDialog()
     <html>
         <head/>
         <body>
-            <p>This is <b>Cell</b>, version 1.1.0. Build time: %1 %2.</p>
+            <p>This is <b>Cell</b>, version 1.1.1. Build time: %1 %2.</p>
             <p>Application for simulating reaction networks.</p>
             <p>Developed by: Yannik Schroeder</p>
             <p>For additional information, see: <a href='https://github.com/yannik131/Cell'>https://github.com/yannik131/Cell</a></p>

src/lib/gui/core/PlotCategories.cpp

@@ -12,7 +12,9 @@ const std::vector<Pair> PlotCategoryPairs{
     {PlotCategory::CollisionCounts, "Collision counts"},
     {PlotCategory::AbsoluteMomentum, "Absolute momentum"},
     {PlotCategory::KineticEnergy, "Kinetic energy"},
-    {PlotCategory::VelocityDistribution, "Velocity distribution"},
+    {PlotCategory::XVelocityDistribution, "X Velocity distribution"},
+    {PlotCategory::YVelocityDistribution, "Y Velocity distribution"},
+    {PlotCategory::AbsoluteVelocityDistribution, "|Velocity| distribution"},
     {PlotCategory::VelocityColorMap, "Velocity color map"},
 };
 

src/lib/gui/core/PlotCategories.hpp

@@ -13,7 +13,9 @@ enum class PlotCategory
     CollisionCounts,
     AbsoluteMomentum,
     KineticEnergy,
-    VelocityDistribution,
+    XVelocityDistribution,
+    YVelocityDistribution,
+    AbsoluteVelocityDistribution,
     VelocityColorMap
 };
 

src/lib/gui/core/Simulation.cpp

@@ -54,7 +54,7 @@ void Simulation::emitLastFrame()
     if (!simulationRecorder_)
         return;
 
-    emit initialFrame(simulationRecorder_->getLastFrame());
+    emit drawFrameImmediately(simulationRecorder_->getLastFrame());
 }
 
 SimulationConfigUpdater& Simulation::getSimulationConfigUpdater()

src/lib/gui/core/Simulation.hpp

@@ -45,7 +45,7 @@ class Simulation : public QObject
 signals:
     void started();
     void stopped();
-    void initialFrame(Frame frame);
+    void drawFrameImmediately(Frame frame);
     void frame(Frame frame);
     void performanceData(const cell::SimulationRunner::PerformanceData& performanceData);
     void dataPoint(const cell::DataPoint& dataPoint);

src/lib/gui/dialogs/SetupDialog.cpp

@@ -3,6 +3,8 @@
 #include "models/SetupModel.hpp"
 #include "ui_SetupDialog.h"
 
+#include <QCheckBox>
+
 SetupDialog::SetupDialog(QWidget* parent, SimulationConfigUpdater* simulationConfigUpdater)
     : QDialog(parent)
     , ui(new Ui::SetupDialog)
@@ -21,6 +23,8 @@ SetupDialog::SetupDialog(QWidget* parent, SimulationConfigUpdater* simulationCon
     connect(ui->cancelPushButton, &QPushButton::clicked, this, &QDialog::reject);
 
     connect(ui->useDistributionRadioButton, &QRadioButton::toggled, setupModel_, &SetupModel::setUseDistribution);
+    connect(ui->conserveAreaCheckbox, &QCheckBox::checkStateChanged, setupModel_,
+            &SetupModel::setReactionsConserveArea);
     connect(ui->timeStepSpinBox, &QSpinBox::valueChanged, setupModel_, &SetupModel::setTimeStepUs);
     connect(ui->timeScaleDoubleSpinBox, &QDoubleSpinBox::valueChanged, setupModel_, &SetupModel::setTimeScale);
     connect(ui->mostProbableSpeedSpinBox, &QDoubleSpinBox::valueChanged, setupModel_,
@@ -48,6 +52,7 @@ void SetupDialog::displayCurrentConfig()
 {
     const auto& config = simulationConfigUpdater_->getSimulationConfig();
     ui->useDistributionRadioButton->setChecked(config.useDistribution);
+    ui->conserveAreaCheckbox->setChecked(config.reactionsConserveArea);
     ui->manualPositionsRadioButton->setChecked(!config.useDistribution);
     ui->timeStepSpinBox->setValue(static_cast<int>(std::round(config.simulationTimeStep * 1e6)));
     ui->timeScaleDoubleSpinBox->setValue(config.simulationTimeScale);

src/lib/gui/dialogs/SetupDialog.ui

@@ -10,7 +10,7 @@
     <x>0</x>
     <y>0</y>
     <width>517</width>
-    <height>325</height>
+    <height>331</height>
    </rect>
   </property>
   <property name="windowTitle">
@@ -183,6 +183,26 @@
         </layout>
        </widget>
       </item>
+      <item>
+       <widget class="QWidget" name="widget_7" native="true">
+        <layout class="QHBoxLayout" name="horizontalLayout_7">
+         <item>
+          <widget class="QLabel" name="label_4">
+           <property name="text">
+            <string>Conserve area in reactions</string>
+           </property>
+          </widget>
+         </item>
+         <item>
+          <widget class="QCheckBox" name="conserveAreaCheckbox">
+           <property name="text">
+            <string/>
+           </property>
+          </widget>
+         </item>
+        </layout>
+       </widget>
+      </item>
      </layout>
     </widget>
    </item>