Add interface file

.gitignore

@@ -123,5 +123,8 @@ doc/latex
 # Draft code for comparing models
 draft_code_my_own_analysis/
 
+# Python environment
+.venv/
+
 # Temporary files
 cout.txt

params.json

@@ -0,0 +1,10 @@
+{
+  "seed": 1,
+  "output_performance_statistics": true,
+  "total_minutes_to_simulate": 43200,
+  "initial_tumor_radius": 150.0,
+  "treatment": {
+    "0": 3957,
+    "8": 3957
+  }
+}

src/cart_tumor.cc

@@ -36,28 +36,40 @@
 #include "core/resource_manager.h"
 #include "core/scheduler.h"
 #include "core/simulation.h"
+#include <cstdint>
 #include <iostream>
 #include <memory>
 #include <vector>
 
 namespace bdm {
 
 int Simulate(int argc, const char** argv) {
-  // Set simulation bounds
-  auto set_param = [](Param* param) {
-    // Set a fixed random seed for reproducibility
-    param->random_seed = kSeed;
-    // Periodic boundary conditions
+  // Load parameters from JSON file or use default values
+  std::unique_ptr<SimParam> custom_parameters = std::make_unique<SimParam>();
+  custom_parameters->LoadParams("params.json");
+
+  // Keep a reference to the parameters before releasing
+  const auto* sparam_ref = custom_parameters.get();
+
+  // Transfer ownership to BioDynaMo parameter system
+  Param::RegisterParamGroup(custom_parameters.release());
+
+  // Set simulation parameters using lambda
+  auto set_param = [sparam_ref](Param* param) {
+    // Set simulation bounds using the parameters
+    param->random_seed = sparam_ref->seed;
     param->bound_space = Param::BoundSpaceMode::kTorus;
-    // Cube of kBoundedSpaceLength³ centered at origin
-    param->min_bound = -kBoundedSpaceLength / kHalf;
-    param->max_bound = kBoundedSpaceLength / kHalf;
-    param->simulation_time_step = kDt;
-    // for outputing performance statistics
-    param->statistics = kOutputPerformanceStatistics;
+    param->min_bound = -sparam_ref->bounded_space_length / kHalf;
+    param->max_bound = sparam_ref->bounded_space_length / kHalf;
+    param->simulation_time_step = sparam_ref->dt_step;
+    param->statistics = sparam_ref->output_performance_statistics;
   };
 
   Simulation simulation(argc, argv, set_param);
+  const auto* sparam = simulation.GetParam()->Get<SimParam>();
+  // Print parameters
+  sparam->PrintParams();
+
   ExecutionContext* ctxt = simulation.GetExecutionContext();
 
   // Change Forces
@@ -72,7 +84,7 @@ int Simulate(int argc, const char** argv) {
   auto* env = dynamic_cast<UniformGridEnvironment*>(
       Simulation::GetActive()->GetEnvironment());
   // Fix the box length for the uniform grid environment
-  env->SetBoxLength(gKLengthBoxMechanics);
+  env->SetBoxLength(sparam->length_box_mechanics);
 
   // Define Substances
   ResourceManager* rm = Simulation::GetActive()->GetResourceManager();
@@ -82,8 +94,9 @@ int Simulate(int argc, const char** argv) {
   // time_step
   std::unique_ptr<DiffusionThomasAlgorithm> oxygen_grid =
       std::make_unique<DiffusionThomasAlgorithm>(
-          kOxygen, "oxygen", kDiffusionCoefficientOxygen, kDecayConstantOxygen,
-          kResolutionGridSubstances, kDtSubstances,
+          kOxygen, "oxygen", sparam->diffusion_coefficient_oxygen,
+          sparam->decay_constant_oxygen, sparam->resolution_grid_substances,
+          sparam->dt_substances,
           /*dirichlet_border=*/true);
   rm->AddContinuum(oxygen_grid.release());
 
@@ -92,9 +105,9 @@ int Simulate(int argc, const char** argv) {
   std::unique_ptr<DiffusionThomasAlgorithm> immunostimulatory_factor_grid =
       std::make_unique<DiffusionThomasAlgorithm>(
           kImmunostimulatoryFactor, "immunostimulatory_factor",
-          kDiffusionCoefficientImmunostimulatoryFactor,
-          kDecayConstantImmunostimulatoryFactor, kResolutionGridSubstances,
-          kDtSubstances,
+          sparam->diffusion_coefficient_immunostimulatory_factor,
+          sparam->decay_constant_immunostimulatory_factor,
+          sparam->resolution_grid_substances, sparam->dt_substances,
           /*dirichlet_border=*/false);
   rm->AddContinuum(immunostimulatory_factor_grid.release());
 
@@ -103,9 +116,10 @@ int Simulate(int argc, const char** argv) {
   // Oxygen comming from the borders (capillary vessels)
   ModelInitializer::AddBoundaryConditions(
       kOxygen, BoundaryConditionType::kDirichlet,
-      // kOxygenReferenceLevel mmHg is the physiological level of oxygen in
+      // oxygen_reference_level mmHg is the physiological level of oxygen in
       // tissues, o2 saturation is 100% at this level
-      std::make_unique<ConstantBoundaryCondition>(kOxygenReferenceLevel));
+      std::make_unique<ConstantBoundaryCondition>(
+          sparam->oxygen_reference_level));
 
   // This is useless now but should be added this way in a future version of
   // BioDynaMo
@@ -114,15 +128,15 @@ int Simulate(int argc, const char** argv) {
 
   // Initialize oxygen voxels
   ModelInitializer::InitializeSubstance(
-      kOxygen, [](real_t /*x*/, real_t /*y*/, real_t /*z*/) {
-        // Set all voxels to kInitialOxygenLevel mmHg
-        return kInitialOxygenLevel;
+      kOxygen, [sparam](real_t /*x*/, real_t /*y*/, real_t /*z*/) {
+        // Set all voxels to initial_oxygen_level mmHg
+        return sparam->initial_oxygen_level;
       });
 
-  // One spherical tumor of radius kInitialRadiusTumor in the center of the
+  // One spherical tumor of radius initial_tumor_radius in the center of the
   // simulation space
   const std::vector<Real3> positions =
-      CreateSphereOfTumorCells(kInitialRadiusTumor);
+      CreateSphereOfTumorCells(sparam->initial_tumor_radius);
   for (const auto& pos : positions) {
     std::unique_ptr<TumorCell> tumor_cell = std::make_unique<TumorCell>(pos);
     std::unique_ptr<StateControlGrowProliferate> state_control =
@@ -133,23 +147,26 @@ int Simulate(int argc, const char** argv) {
 
   // Treatment administration operation
   std::unique_ptr<bdm::Operation> treatment_op =
-      std::make_unique<bdm::Operation>("SpawnCart", kStepsOneDay);
+      std::make_unique<bdm::Operation>("SpawnCart", sparam->steps_in_one_day);
   std::unique_ptr<bdm::SpawnCart> spawn_cart =
       std::make_unique<bdm::SpawnCart>();
   treatment_op->AddOperationImpl(bdm::kCpu, spawn_cart.release());
   scheduler->ScheduleOp(treatment_op.release());
 
   // OutputSummary operation
-  std::unique_ptr<bdm::Operation> summary_op =
-      std::make_unique<bdm::Operation>("OutputSummary", kOutputCsvInterval);
+  std::unique_ptr<bdm::Operation> summary_op = std::make_unique<bdm::Operation>(
+      "OutputSummary", sparam->output_csv_interval);
   std::unique_ptr<bdm::OutputSummary> output_summary =
       std::make_unique<bdm::OutputSummary>();
   summary_op->AddOperationImpl(bdm::kCpu, output_summary.release());
   scheduler->ScheduleOp(summary_op.release());
 
   // Run simulation
-  // simulate kTotalMinutesToSimulate minutes including the last minute
-  scheduler->Simulate(1 + kTotalMinutesToSimulate / kDt);
+  std::cout << "Running simulation..." << std::endl;
+  // simulate total_minutes_to_simulate minutes including the last minute
+  scheduler->Simulate(1 +
+                      static_cast<uint64_t>(sparam->total_minutes_to_simulate /
+                                            sparam->dt_step));
   std::cout << "Simulation completed successfully!" << std::endl;
   return 0;
 }

src/diffusion_thomas_algorithm.cc

@@ -26,6 +26,7 @@
 #include "core/agent/agent.h"
 #include "core/container/math_array.h"
 #include "core/diffusion/diffusion_grid.h"
+#include "core/param/param.h"
 #include "core/real_t.h"
 #include "core/resource_manager.h"
 #include <cstddef>
@@ -44,12 +45,15 @@ DiffusionThomasAlgorithm::DiffusionThomasAlgorithm(int substance_id,
     : DiffusionGrid(substance_id, std::move(substance_name), dc, mu,
                     resolution),
       resolution_(static_cast<int>(GetResolution())),
-      d_space_(static_cast<real_t>(kBoundedSpaceLength) /
+      d_space_(static_cast<real_t>(Simulation::GetActive()
+                                       ->GetParam()
+                                       ->Get<SimParam>()
+                                       ->bounded_space_length) /
                static_cast<real_t>(resolution_)),
       dirichlet_border_(dirichlet_border),
       jump_i_(1),
       jump_j_(resolution_),
-      jump_k_(resolution_ * resolution_),
+      jump_(resolution_ * resolution_),
       spatial_diffusion_coeff_(dc * dt / (d_space_ * d_space_)),
       neg_diffusion_factor_(-spatial_diffusion_coeff_),
       temporal_decay_coeff_(mu * dt / 3.0),
@@ -63,7 +67,6 @@ DiffusionThomasAlgorithm::DiffusionThomasAlgorithm(int substance_id,
       thomas_c_z_(resolution_, neg_diffusion_factor_),
       thomas_denom_z_(resolution_, central_coeff_) {
   SetTimeStep(dt);
-
   // Initialize the denominators and coefficients for the Thomas algorithm
   InitializeThomasAlgorithmVectors(thomas_denom_x_, thomas_c_x_);
   InitializeThomasAlgorithmVectors(thomas_denom_y_, thomas_c_y_);
@@ -245,7 +248,7 @@ void DiffusionThomasAlgorithm::SolveDirectionThomas(int direction) {
       return jump_j_;
     }
     // direction == 2
-    return jump_k_;
+    return jump_;
   }();
 
 #pragma omp parallel for collapse(2)

src/diffusion_thomas_algorithm.h

@@ -43,7 +43,7 @@ class DiffusionThomasAlgorithm : public DiffusionGrid {
         dirichlet_border_(false),
         jump_i_(0),
         jump_j_(0),
-        jump_k_(0),
+        jump_(0),
         spatial_diffusion_coeff_(0.0),
         neg_diffusion_factor_(0.0),
         temporal_decay_coeff_(0.0),
@@ -138,7 +138,7 @@ class DiffusionThomasAlgorithm : public DiffusionGrid {
   int jump_j_;
 
   /// Index jump for k-direction (z-axis)
-  int jump_k_;
+  int jump_;
 
   /// First diffusion constant
   real_t spatial_diffusion_coeff_;

src/forces_tumor_cart.cc

@@ -27,6 +27,7 @@
 #include "core/agent/cell.h"
 #include "core/container/math_array.h"
 #include "core/interaction_force.h"
+#include "core/param/param.h"
 #include "core/real_t.h"
 #include <algorithm>
 #include <cmath>
@@ -43,18 +44,23 @@ Real4 InteractionVelocity::Calculate(const Agent* lhs, const Agent* rhs) const {
     return {0.0, 0.0, 0.0, 0.0};
   }
 
+  const auto* sparams = Simulation::GetActive()->GetParam()->Get<SimParam>();
+
   Real3 displacement = a->GetPosition() - b->GetPosition();
 
   // For periodic boundary conditions, we need to adjust the displacement
   displacement[0] =
       displacement[0] -
-      (kBoundedSpaceLength)*round(displacement[0] / (kBoundedSpaceLength));
+      (sparams->bounded_space_length) *
+          round(displacement[0] / (sparams->bounded_space_length));
   displacement[1] =
       displacement[1] -
-      (kBoundedSpaceLength)*round(displacement[1] / (kBoundedSpaceLength));
+      (sparams->bounded_space_length) *
+          round(displacement[1] / (sparams->bounded_space_length));
   displacement[2] =
       displacement[2] -
-      (kBoundedSpaceLength)*round(displacement[2] / (kBoundedSpaceLength));
+      (sparams->bounded_space_length) *
+          round(displacement[2] / (sparams->bounded_space_length));
 
   const real_t dist_sq = displacement[0] * displacement[0] +
                          displacement[1] * displacement[1] +
@@ -79,23 +85,25 @@ Real4 InteractionVelocity::Calculate(const Agent* lhs, const Agent* rhs) const {
 
     if ((a_tumor != nullptr) && (b_tumor != nullptr)) {
       // two tumor cells
-      // std::sqrt(kRepulsionTumorTumor * kRepulsionTumorTumor);
-      repulsion = kRepulsionTumorTumor;
+      // std::sqrt(sparams->cell_repulsion_between_tumor_tumor *
+      // sparams->cell_repulsion_between_tumor_tumor);
+      repulsion = sparams->cell_repulsion_between_tumor_tumor;
     } else if ((a_tumor == nullptr) && (b_tumor == nullptr)) {
       // two CAR-T cells
-      // std::sqrt(kRepulsionCartCart*kRepulsionCartCart);
-      repulsion = kRepulsionCartCart;
+      // std::sqrt(sparams->cell_repulsion_between_cart_cart*sparams->cell_repulsion_between_cart_cart);
+      repulsion = sparams->cell_repulsion_between_cart_cart;
     } else {
       // one tumor cell and one CAR-T
-      repulsion = std::sqrt(kRepulsionCartTumor * kRepulsionTumorCart);
+      repulsion = std::sqrt(sparams->cell_repulsion_between_cart_tumor *
+                            sparams->cell_repulsion_between_tumor_cart);
     }
 
     temp_r *= repulsion;
   }
 
   // Adhesion
   const real_t max_interaction_distance =
-      kMaxRelativeAdhesionDistance * combined_radius;
+      sparams->max_relative_adhesion_distance * combined_radius;
 
   if (distance < max_interaction_distance) {
     // 1 - d/S
@@ -106,13 +114,14 @@ Real4 InteractionVelocity::Calculate(const Agent* lhs, const Agent* rhs) const {
     real_t adhesion = NAN;
     if ((a_tumor != nullptr) && (b_tumor != nullptr)) {
       // two tumor cells
-      adhesion = kAdhesionTumorTumor;
+      adhesion = sparams->cell_adhesion_between_tumor_tumor;
     } else if ((a_tumor == nullptr) && (b_tumor == nullptr)) {
       // two CAR-T cells
-      adhesion = kAdhesionCartCart;
+      adhesion = sparams->cell_adhesion_between_cart_cart;
     } else {
       // one tumor cell and one CAR-T
-      adhesion = std::sqrt(kAdhesionCartTumor * kAdhesionTumorCart);
+      adhesion = std::sqrt(sparams->cell_adhesion_between_cart_tumor *
+                           sparams->cell_adhesion_between_tumor_cart);
     }
 
     temp_a *= adhesion;