diff --git a/demo/analytic_continuum/src/continuum_model.h b/demo/analytic_continuum/src/continuum_model.h index 1ea5ab146..f984e0866 100644 --- a/demo/analytic_continuum/src/continuum_model.h +++ b/demo/analytic_continuum/src/continuum_model.h @@ -26,7 +26,6 @@ namespace bdm { class AnalyticContinuum final : public ScalarField { public: AnalyticContinuum() = default; - explicit AnalyticContinuum(const TRootIOCtor *) {} ~AnalyticContinuum() final = default; void Initialize() final {} @@ -58,8 +57,6 @@ class AnalyticContinuum final : public ScalarField { return {0, 0, 0}; }; - BDM_CLASS_DEF_OVERRIDE(AnalyticContinuum, 1); // NOLINT - private: real_t time_ = 0.0; }; diff --git a/demo/analytic_continuum/src/my_agent.h b/demo/analytic_continuum/src/my_agent.h index 23e42d89c..f459e8177 100644 --- a/demo/analytic_continuum/src/my_agent.h +++ b/demo/analytic_continuum/src/my_agent.h @@ -23,7 +23,7 @@ namespace bdm { /// A simple agent that can sense the continuum value and stores the value in /// a member variable. class ContinuumRetrieverAgent : public SphericalAgent { - BDM_AGENT_HEADER(ContinuumRetrieverAgent, SphericalAgent, 1); + BDM_AGENT_HEADER(ContinuumRetrieverAgent, SphericalAgent); public: ContinuumRetrieverAgent() {} diff --git a/demo/analytic_continuum/src/my_behavior.h b/demo/analytic_continuum/src/my_behavior.h index 89771bba0..7229cb72e 100644 --- a/demo/analytic_continuum/src/my_behavior.h +++ b/demo/analytic_continuum/src/my_behavior.h @@ -24,7 +24,7 @@ namespace bdm { /// Behavior that allows agents to sense the continuum value and stores the /// value in its member variable. struct RetrieveContinuumValue : public Behavior { - BDM_BEHAVIOR_HEADER(RetrieveContinuumValue, Behavior, 1); + BDM_BEHAVIOR_HEADER(RetrieveContinuumValue, Behavior); RetrieveContinuumValue() {} virtual ~RetrieveContinuumValue() {} diff --git a/demo/cell_division_enhanced/src/cell_division_enhanced.h b/demo/cell_division_enhanced/src/cell_division_enhanced.h index 21205633d..76e944a3f 100644 --- a/demo/cell_division_enhanced/src/cell_division_enhanced.h +++ b/demo/cell_division_enhanced/src/cell_division_enhanced.h @@ -33,7 +33,7 @@ inline int signum(const T& x) { namespace bdm { class MyCell : public Cell { - BDM_AGENT_HEADER(MyCell, Cell, 1); + BDM_AGENT_HEADER(MyCell, Cell); public: MyCell() : Cell() { UpdateVolume(); } diff --git a/demo/flocking/src/boid.h b/demo/flocking/src/boid.h index 2d5dadf36..385867107 100644 --- a/demo/flocking/src/boid.h +++ b/demo/flocking/src/boid.h @@ -50,7 +50,7 @@ real_t FluctuateCoefficient(real_t coefficient, real_t fluctuation_strength, //////////////////////////////////////////////////////////////////////////////// class Boid : public Agent { - BDM_AGENT_HEADER(Boid, Agent, 1); + BDM_AGENT_HEADER(Boid, Agent); public: Boid() {} @@ -171,13 +171,13 @@ class Boid : public Agent { //////////////////////////////////////////////////////////////////////////////// struct Flocking : public Behavior { - BDM_BEHAVIOR_HEADER(Flocking, Behavior, 1); + BDM_BEHAVIOR_HEADER(Flocking, Behavior); void Run(Agent* agent) override; }; struct RandomPerturbation : public Behavior { - BDM_BEHAVIOR_HEADER(RandomPerturbation, Behavior, 1); + BDM_BEHAVIOR_HEADER(RandomPerturbation, Behavior); public: RandomPerturbation(real_t v = 1) : velocity_(v) { AlwaysCopyToNew(); } diff --git a/demo/mars/src/CelestialObjects.h b/demo/mars/src/CelestialObjects.h index 72a71d41a..c6ca29859 100644 --- a/demo/mars/src/CelestialObjects.h +++ b/demo/mars/src/CelestialObjects.h @@ -22,7 +22,7 @@ namespace astrophysics { // creating the custom agent that // defines celestial objects class CelestialObject : public Cell { - BDM_AGENT_HEADER(CelestialObject, Cell, 1); + BDM_AGENT_HEADER(CelestialObject, Cell); public: // constructors @@ -51,7 +51,7 @@ class CelestialObject : public Cell { // Planet subclass class Planet : public CelestialObject { - BDM_AGENT_HEADER(Planet, CelestialObject, 1); + BDM_AGENT_HEADER(Planet, CelestialObject); public: Planet() : CelestialObject() {} @@ -65,7 +65,7 @@ class Planet : public CelestialObject { // Satellite subclass class Satellite : public CelestialObject { - BDM_AGENT_HEADER(Satellite, CelestialObject, 1); + BDM_AGENT_HEADER(Satellite, CelestialObject); public: Satellite() : CelestialObject() {} diff --git a/demo/monolayer_growth/src/behaviours.h b/demo/monolayer_growth/src/behaviours.h index aeb94632d..df5d887a8 100644 --- a/demo/monolayer_growth/src/behaviours.h +++ b/demo/monolayer_growth/src/behaviours.h @@ -24,7 +24,7 @@ namespace bdm { // Define growth behaviour struct GrowthAndCellCycle : public Behavior { - BDM_BEHAVIOR_HEADER(GrowthAndCellCycle, Behavior, 1); + BDM_BEHAVIOR_HEADER(GrowthAndCellCycle, Behavior); GrowthAndCellCycle() { AlwaysCopyToNew(); } virtual ~GrowthAndCellCycle() {} diff --git a/demo/monolayer_growth/src/cell_cell_force.cc b/demo/monolayer_growth/src/cell_cell_force.cc index 90972857f..ef01cf835 100644 --- a/demo/monolayer_growth/src/cell_cell_force.cc +++ b/demo/monolayer_growth/src/cell_cell_force.cc @@ -17,6 +17,10 @@ namespace bdm { +namespace { +constexpr real_t kCoincidentCenterTolerance = 1e-8; +} // namespace + /// Custom force. Changed adhesive and repulsive parameters compared to standard /// force to achieve quick separation of mother and daughter cells after /// division. @@ -52,7 +56,7 @@ Real4 CellCellForce::Calculate(const Agent* lhs, const Agent* rhs) const { } // to avoid a division by 0 if the centers are (almost) at the same // location - if (center_distance < 0.00000001) { + if (center_distance < kCoincidentCenterTolerance) { auto* random = Simulation::GetActive()->GetRandom(); auto force2on1 = random->template UniformArray<3>(-3.0, 3.0); return {force2on1[0], force2on1[1], force2on1[2], 0}; diff --git a/demo/monolayer_growth/src/cycling_cell.h b/demo/monolayer_growth/src/cycling_cell.h index 86dda0361..4cbb34ef9 100644 --- a/demo/monolayer_growth/src/cycling_cell.h +++ b/demo/monolayer_growth/src/cycling_cell.h @@ -25,7 +25,7 @@ enum CellState { kG1, kS, kG2, kM }; class CyclingCell : public Cell { // our object extends the Cell object // create the header with our new data member - BDM_AGENT_HEADER(CyclingCell, Cell, 1); + BDM_AGENT_HEADER(CyclingCell, Cell); public: CyclingCell() {} diff --git a/demo/monolayer_growth/src/evaluate.h b/demo/monolayer_growth/src/evaluate.h index ce46ed72a..8a4469cbb 100644 --- a/demo/monolayer_growth/src/evaluate.h +++ b/demo/monolayer_growth/src/evaluate.h @@ -15,16 +15,13 @@ #ifndef EVALUATE_H_ #define EVALUATE_H_ -#include -#include -#include #include +#include +#include #include #include "biodynamo.h" #include "sim_param.h" -using namespace bdm::experimental; - namespace bdm { inline void SetupResultCollection(Simulation* sim) { @@ -56,46 +53,27 @@ inline void SetupResultCollection(Simulation* sim) { ts->AddCollector("env_dims", get_env_dims, get_time); } -inline void ExportResults(const bool plot_legend = true, - const std::string& filename = "result") { - // Prerequisites +inline void ExportResults(const std::string& filename = "result.csv") { const std::string folder = Simulation::GetActive()->GetOutputDir(); auto* ts = Simulation::GetActive()->GetTimeSeries(); - TimeSeries allts; - std::vector times, sizes; - - // Add simulated data - allts.Add(*ts, Concat("i", 0)); - times = ts->GetXValues("env_dims"); - sizes = ts->GetYValues("env_dims"); - - // Add experimental data from Figure 1 from Drasdo and Hoehme (2005) - allts.Add("experimental_data", - {336 / 24., 386 / 24., 408 / 24., 481 / 24., 506 / 24., 646 / 24.}, - {1140, 1400, 1590, 2040, 2250, 3040}); - - // Initialize line graph - LineGraph lg(&allts, "", "Time [days]", "2D Monolayer size (um)", plot_legend, - nullptr, 350, 250); - - // Add simulated data - lg.Add(Concat("env_dims-i", 0), "Sim data ", "LP", kBlue, 0.2, kSolid, 2, - kBlue, 0.7, kFullCircle, 0.5); - - // Add style for exp data - lg.Add("experimental_data", "Exp data ", "LP", kBlack, 0.2, kSolid, 2, kBlack, - 0.7, kFullCircle, 0.5); - - // Add legend - if (plot_legend) { - lg.SetLegendPosNDC(0.1, 0.7, 0.3, 0.9); + const auto& times = ts->GetXValues("env_dims"); + const auto& sizes = ts->GetYValues("env_dims"); + std::ofstream output(Concat(folder, "/", filename)); + output << "source,time_days,size_um\n"; + for (size_t i = 0; i < times.size(); ++i) { + output << "simulation," << times[i] << ',' << sizes[i] << '\n'; } - // Save plot - lg.SaveAs(Concat(folder, "/", filename), {".svg", ".png"}); - - // Save data - ts->SaveJson(Concat(folder, "/monolayer_growth.json")); + constexpr real_t kHoursPerDay = 24.0; + const std::vector experimental_times = { + 336 / kHoursPerDay, 386 / kHoursPerDay, 408 / kHoursPerDay, + 481 / kHoursPerDay, 506 / kHoursPerDay, 646 / kHoursPerDay}; + const std::vector experimental_sizes = {1140, 1400, 1590, + 2040, 2250, 3040}; + for (size_t i = 0; i < experimental_times.size(); ++i) { + output << "experiment," << experimental_times[i] << ',' + << experimental_sizes[i] << '\n'; + } } } // namespace bdm diff --git a/demo/newtons_law_test/src/CelestialObjects.h b/demo/newtons_law_test/src/CelestialObjects.h index 7e148cdcc..66e892f13 100644 --- a/demo/newtons_law_test/src/CelestialObjects.h +++ b/demo/newtons_law_test/src/CelestialObjects.h @@ -22,7 +22,7 @@ namespace astrophysics { // creating the custom agent that // defines celestial objects class CelestialObject : public Cell { - BDM_AGENT_HEADER(CelestialObject, Cell, 1); + BDM_AGENT_HEADER(CelestialObject, Cell); public: // constructors diff --git a/demo/parameters/src/parameters.h b/demo/parameters/src/parameters.h index ba5ce8752..1ef485dbb 100644 --- a/demo/parameters/src/parameters.h +++ b/demo/parameters/src/parameters.h @@ -24,6 +24,13 @@ struct SimParam : public ParamGroup { real_t foo = 3.14; int bar = -42; + + protected: + void AssignFromConfig( + const std::shared_ptr& config) override { + BDM_ASSIGN_CONFIG_VALUE(foo, "parameters.foo"); + BDM_ASSIGN_CONFIG_VALUE(bar, "parameters.bar"); + } }; inline int Simulate(int argc, const char** argv) { diff --git a/demo/pyramidal_cell/src/pyramidal_cell.h b/demo/pyramidal_cell/src/pyramidal_cell.h index 8ff85b740..d644cd532 100644 --- a/demo/pyramidal_cell/src/pyramidal_cell.h +++ b/demo/pyramidal_cell/src/pyramidal_cell.h @@ -24,7 +24,7 @@ namespace bdm { enum Substances { kApical, kBasal }; struct ApicalDendriteGrowth : public Behavior { - BDM_BEHAVIOR_HEADER(ApicalDendriteGrowth, Behavior, 1); + BDM_BEHAVIOR_HEADER(ApicalDendriteGrowth, Behavior); ApicalDendriteGrowth() { AlwaysCopyToNew(); } virtual ~ApicalDendriteGrowth() {} @@ -83,7 +83,7 @@ struct ApicalDendriteGrowth : public Behavior { }; struct BasalDendriteGrowth : public Behavior { - BDM_BEHAVIOR_HEADER(BasalDendriteGrowth, Behavior, 1); + BDM_BEHAVIOR_HEADER(BasalDendriteGrowth, Behavior); BasalDendriteGrowth() { AlwaysCopyToNew(); } virtual ~BasalDendriteGrowth() {} diff --git a/demo/regulatory_networks/src/bdm_ex1.h b/demo/regulatory_networks/src/bdm_ex1.h index 780f9f150..fab00f8d7 100644 --- a/demo/regulatory_networks/src/bdm_ex1.h +++ b/demo/regulatory_networks/src/bdm_ex1.h @@ -22,7 +22,6 @@ namespace bdm { enum Substances { kProtein }; -#ifndef __ROOTCLING__ struct ODE_system { const std::map& mdg; const std::vector param; @@ -80,8 +79,6 @@ struct ODE_output { std::clog << std::endl; } }; -#endif - namespace ex1 { inline int Simulate(int argc, const char** argv) { @@ -136,15 +133,12 @@ inline int Simulate(int argc, const char** argv) { c->SetAdherence(0.4); c->SetMass(1.0); c->SetPosition(xyz); -#ifndef __ROOTCLING__ c->AddBehavior(new RegulatoryNetwork( dt_RN, 1000, {1., 5., 7.}, // ODE_solver::Euler, // ODE_solver::Rosenbrock, ODE_solver::RungeKutta, ODE_system(dg_map, {0.2, 0.1, 3.0}), ODE_jacobian(dg_map, {0.2, 0.1, 3.0}), ODE_output())); -#endif - sim.GetExecutionContext()->AddAgent(c); } diff --git a/demo/regulatory_networks/src/bdm_ex2.h b/demo/regulatory_networks/src/bdm_ex2.h index 610d0adfb..88a7999d6 100644 --- a/demo/regulatory_networks/src/bdm_ex2.h +++ b/demo/regulatory_networks/src/bdm_ex2.h @@ -22,7 +22,7 @@ namespace bdm { class MyCell : public Cell { - BDM_AGENT_HEADER(MyCell, Cell, 1); + BDM_AGENT_HEADER(MyCell, Cell); public: MyCell() {} @@ -42,12 +42,20 @@ class MyCell : public Cell { real_t GetTrail() const { return trail_; } void SetTrail(real_t t) { trail_ += t; } + bool GetVisualizationData(const std::string& name, + VisualizationData* values) const override { + if (name == "trail_") { + *values = std::vector{trail_}; + return true; + } + return Base::GetVisualizationData(name, values); + } + private: /// keep track of the trail of the agent real_t trail_; }; -#ifndef __ROOTCLING__ struct Lorenz_rhs_ { void operator()(const b_vector_t& x, b_vector_t& dxdt, double t, Agent* agent) const { @@ -91,18 +99,14 @@ struct Lorenz_out_ { std::clog << std::endl; } }; -#endif - class Trajectory : public RegulatoryNetwork { - BDM_BEHAVIOR_HEADER(Trajectory, RegulatoryNetwork, 1); + BDM_BEHAVIOR_HEADER(Trajectory, RegulatoryNetwork); public: Trajectory() { AlwaysCopyToNew(); } -#ifndef __ROOTCLING__ Trajectory(real_t dt, int n_dt, const std::vector& x) : RegulatoryNetwork(dt, n_dt, x, ODE_solver::Rosenbrock, Lorenz_rhs_(), Lorenz_jac_(), Lorenz_out_()) {} -#endif virtual ~Trajectory() = default; void Initialize(const NewAgentEvent& event) override { @@ -112,7 +116,6 @@ class Trajectory : public RegulatoryNetwork { void Run(Agent* agent) override { Base::Run(agent); -#ifndef __ROOTCLING__ Real3 xyz; for (int i = 0; i < 3; i++) xyz[i] = this->GetSpecie(i); @@ -126,7 +129,6 @@ class Trajectory : public RegulatoryNetwork { } else { Log::Fatal("Trajectory::Run", "agent is not of 'MyCell' type"); } -#endif } }; @@ -164,11 +166,7 @@ inline int Simulate(int argc, const char** argv) { MyCell* c = new MyCell(); c->SetDiameter(1.0); c->SetPosition(xyz); -#ifndef __ROOTCLING__ c->AddBehavior(new Trajectory(dt_RN, 222, {xyz[0], xyz[1], xyz[2]})); -#else - c->AddBehavior(new Trajectory()); -#endif sim.GetExecutionContext()->AddAgent(c); } diff --git a/demo/regulatory_networks/src/lorenz.h b/demo/regulatory_networks/src/lorenz.h index 309ffc274..685710e3d 100644 --- a/demo/regulatory_networks/src/lorenz.h +++ b/demo/regulatory_networks/src/lorenz.h @@ -18,21 +18,16 @@ #include #include -#ifndef __ROOTCLING__ #include #include "boost/numeric/odeint.hpp" #include "boost/phoenix/core.hpp" #include "boost/phoenix/operator.hpp" -#endif -#ifndef __ROOTCLING__ typedef boost::numeric::ublas::vector b_vector_t; typedef boost::numeric::ublas::matrix b_matrix_t; -#endif namespace lorenz { -#ifndef __ROOTCLING__ struct ODE_system { void operator()(const b_vector_t& x, b_vector_t& dxdt, double t) const { dxdt[0] = sigma * x[1] - sigma * x[0]; @@ -71,8 +66,6 @@ struct ODE_output { std::clog << t << ',' << x[0] << ',' << x[1] << ',' << x[2] << std::endl; } }; -#endif - inline int Simulate(int argc, const char** argv) { std::ofstream fout("lorenz.csv"); // save the original buffer of std::clog @@ -80,7 +73,6 @@ inline int Simulate(int argc, const char** argv) { // redirect std::clog to point to the above file std::clog.rdbuf(fout.rdbuf()); -#ifndef __ROOTCLING__ b_vector_t xyz(3); xyz[0] = 1.0; xyz[1] = 1.0; @@ -94,8 +86,6 @@ inline int Simulate(int argc, const char** argv) { // perform the time-integration integrate_const(stepper, std::make_pair(ODE_system(), ODE_jacobian()), xyz, 0.0, 30.0, 0.01, ODE_output()); -#endif - // restore the original buffer of std::clog std::clog.rdbuf(orig_clog_buff); diff --git a/demo/regulatory_networks/src/oscillator.h b/demo/regulatory_networks/src/oscillator.h index 0dda12b8a..6c02786d3 100644 --- a/demo/regulatory_networks/src/oscillator.h +++ b/demo/regulatory_networks/src/oscillator.h @@ -18,21 +18,16 @@ #include #include -#ifndef __ROOTCLING__ #include #include "boost/numeric/odeint.hpp" #include "boost/phoenix/core.hpp" #include "boost/phoenix/operator.hpp" -#endif -#ifndef __ROOTCLING__ typedef boost::numeric::ublas::vector b_vector_t; typedef boost::numeric::ublas::matrix b_matrix_t; -#endif namespace oscillator { -#ifndef __ROOTCLING__ struct ODE_system { void operator()(const b_vector_t& x, b_vector_t& dxdt, double t) const { dxdt[0] = x[1]; @@ -61,8 +56,6 @@ struct ODE_output { std::clog << t << ',' << x[0] << ',' << x[1] << std::endl; } }; -#endif - inline int Simulate(int argc, const char** argv) { std::ofstream fout("oscillator.csv"); // save the original buffer of std::clog @@ -70,7 +63,6 @@ inline int Simulate(int argc, const char** argv) { // redirect std::clog to point to the above file std::clog.rdbuf(fout.rdbuf()); -#ifndef __ROOTCLING__ b_vector_t xy(2); xy[0] = 2; xy[1] = 0; @@ -83,8 +75,6 @@ inline int Simulate(int argc, const char** argv) { // perform the time-integration integrate_const(stepper, std::make_pair(ODE_system(), ODE_jacobian()), xy, 0.0, 500.0, 1.0e-2, ODE_output()); -#endif - // restore the original buffer of std::clog std::clog.rdbuf(orig_clog_buff); diff --git a/demo/regulatory_networks/src/sine.h b/demo/regulatory_networks/src/sine.h index 3c9aeca01..16a466a14 100644 --- a/demo/regulatory_networks/src/sine.h +++ b/demo/regulatory_networks/src/sine.h @@ -18,21 +18,21 @@ #include #include -#ifndef __ROOTCLING__ #include #include "boost/numeric/odeint.hpp" #include "boost/phoenix/core.hpp" #include "boost/phoenix/operator.hpp" -#endif +#include "core/util/math.h" -#ifndef __ROOTCLING__ typedef boost::numeric::ublas::vector b_vector_t; typedef boost::numeric::ublas::matrix b_matrix_t; -#endif namespace sine { -#ifndef __ROOTCLING__ +constexpr double kIntegrationTolerance = 1e-6; +constexpr double kIntegrationStep = 1e-3; +constexpr double kEndTime = 4.0 * bdm::Math::kPi; + struct ODE_system { void operator()(const b_vector_t& x, b_vector_t& dxdt, double t) const { dxdt[0] = A * cos(t); @@ -46,8 +46,6 @@ struct ODE_output { std::clog << t << ',' << x[0] << std::endl; } }; -#endif - inline int Simulate(int argc, const char** argv) { std::ofstream fout("sine.csv"); // save the original buffer of std::clog @@ -55,20 +53,18 @@ inline int Simulate(int argc, const char** argv) { // redirect std::clog to point to the above file std::clog.rdbuf(fout.rdbuf()); -#ifndef __ROOTCLING__ b_vector_t x(1); x[0] = 0.0; typedef boost::numeric::odeint::runge_kutta_dopri5 ode_int; // set-up the Runge-Kutta integrator - auto stepper = boost::numeric::odeint::make_dense_output(1e-6, 1e-6); + auto stepper = boost::numeric::odeint::make_dense_output( + kIntegrationTolerance, kIntegrationTolerance); // perform the time-integration - integrate_const(stepper, ODE_system(), x, 0.0, 12.5663706144, 0.001, + integrate_const(stepper, ODE_system(), x, 0.0, kEndTime, kIntegrationStep, ODE_output()); -#endif - // restore the original buffer of std::clog std::clog.rdbuf(orig_clog_buff); diff --git a/demo/sbml_integration/src/sbml_integration.h b/demo/sbml_integration/src/sbml_integration.h index f9ef09043..45af1ce72 100644 --- a/demo/sbml_integration/src/sbml_integration.h +++ b/demo/sbml_integration/src/sbml_integration.h @@ -14,17 +14,13 @@ #ifndef SBML_INTEGRATION_H_ #define SBML_INTEGRATION_H_ +#include +#include + #include "biodynamo.h" #include "core/util/io.h" #include "core/util/timing.h" -#include -#include -#include -#include -#include -#include - #include "rrException.h" #include "rrExecutableModel.h" #include "rrLogger.h" @@ -36,7 +32,7 @@ namespace bdm { // Define my custom cell, which extends Cell by adding an extra // data member s1_. class MyCell : public Cell { - BDM_AGENT_HEADER(MyCell, Cell, 1); + BDM_AGENT_HEADER(MyCell, Cell); public: MyCell() {} @@ -52,7 +48,7 @@ class MyCell : public Cell { // Define SbmlBehavior to simulate intracellular chemical reaction network. class SbmlBehavior : public Behavior { - BDM_BEHAVIOR_HEADER(SbmlBehavior, Behavior, 1) + BDM_BEHAVIOR_HEADER(SbmlBehavior, Behavior) public: SbmlBehavior() {} @@ -67,14 +63,14 @@ class SbmlBehavior : public Behavior { result_ = other_sbml_behavior->result_; } - virtual ~SbmlBehavior() { delete rr_; } + virtual ~SbmlBehavior() = default; void Initialize(const std::string& sbml_file, const rr::SimulateOptions& opt) { sbml_file_ = sbml_file; initial_options_ = opt; - rr_ = new rr::RoadRunner(sbml_file); + rr_ = std::make_unique(sbml_file); rr_->getSimulateOptions() = opt; // setup integrator rr_->setIntegrator("gillespie"); @@ -113,66 +109,32 @@ class SbmlBehavior : public Behavior { rr::SimulateOptions initial_options_; ls::DoubleMatrix result_; bool active_ = true; - rr::RoadRunner* rr_; - real_t dt_; + std::unique_ptr rr_; + real_t dt_ = 0; }; -inline void AddToPlot(TMultiGraph* mg, const ls::Matrix* result) { - ls::Matrix foo1(*result); - ls::Matrix foo(*foo1.getTranspose()); - int rows; - int cols; - auto** twod = foo.get2DMatrix(rows, cols); - - TGraph* gr = new TGraph(cols, twod[0], twod[1]); - gr->SetFillStyle(0); - gr->SetLineColorAlpha(2, 0.1); - gr->SetLineWidth(1); - gr->SetTitle("S1"); - - TGraph* gr1 = new TGraph(cols, twod[0], twod[2]); - gr1->SetTitle("S2"); - gr1->SetLineColorAlpha(3, 0.1); - gr1->SetLineWidth(1); - - TGraph* gr2 = new TGraph(cols, twod[0], twod[3]); - gr2->SetTitle("S3"); - gr2->SetLineColorAlpha(4, 0.1); - gr2->SetLineWidth(1); - - mg->Add(gr); - mg->Add(gr1); - mg->Add(gr2); - mg->Draw("AL C C"); -} - -inline void PlotSbmlBehaviors(const char* filename) { - // setup plot - TCanvas c; - c.SetGrid(); - - TMultiGraph* mg = new TMultiGraph(); - mg->SetTitle("Gillespie;Timestep;Concentration"); +inline void ExportSbmlBehaviors(const char* filename) { + std::ofstream output(filename); + output << "agent,time,s1,s2,s3\n"; + uint64_t agent_index = 0; Simulation::GetActive()->GetResourceManager()->ForEachAgent( [&](Agent* agent) { auto* cell = static_cast(agent); const auto& behaviour = cell->GetAllBehaviors(); if (behaviour.size() == 1) { - AddToPlot(mg, &static_cast(behaviour[0])->GetResult()); + const auto& result = + static_cast(behaviour[0])->GetResult(); + for (unsigned row = 0; row < result.numRows(); ++row) { + output << agent_index; + for (unsigned column = 0; column < result.numCols(); ++column) { + output << ',' << result(row, column); + } + output << '\n'; + } } + ++agent_index; }); - - // finalize plot - // TCanvas::Update() draws the frame, after which one can change it - c.Update(); - c.GetFrame()->SetBorderSize(12); - gPad->Modified(); - gPad->Update(); - c.Modified(); - c.cd(0); - // c.BuildLegend(); // TODO position of legend - c.SaveAs(filename); } inline int Simulate(int argc, const char** argv) { @@ -216,7 +178,7 @@ inline int Simulate(int argc, const char** argv) { auto stop = Timing::Timestamp(); std::cout << "RUNTIME " << (stop - start) << std::endl; - PlotSbmlBehaviors("sbml-behaviors.svg"); + ExportSbmlBehaviors("sbml-behaviors.csv"); std::cout << "Simulation completed successfully!" << std::endl; return 0; diff --git a/demo/solar_system/src/CelestialObjects.h b/demo/solar_system/src/CelestialObjects.h index 5bbef099d..36d02cbea 100644 --- a/demo/solar_system/src/CelestialObjects.h +++ b/demo/solar_system/src/CelestialObjects.h @@ -22,7 +22,7 @@ namespace astrophysics { // creating the custom agent that // defines celestial objects class CelestialObject : public Cell { - BDM_AGENT_HEADER(CelestialObject, Cell, 1); + BDM_AGENT_HEADER(CelestialObject, Cell); public: // constructors @@ -51,7 +51,7 @@ class CelestialObject : public Cell { // Star subclass class Star : public CelestialObject { - BDM_AGENT_HEADER(Star, CelestialObject, 1); + BDM_AGENT_HEADER(Star, CelestialObject); public: Star() : CelestialObject() {} @@ -65,7 +65,7 @@ class Star : public CelestialObject { // Planet subclass class Planet : public CelestialObject { - BDM_AGENT_HEADER(Planet, CelestialObject, 1); + BDM_AGENT_HEADER(Planet, CelestialObject); public: Planet() : CelestialObject() {} diff --git a/demo/solar_system/src/solar_system.h b/demo/solar_system/src/solar_system.h index 916c78e95..a81758330 100644 --- a/demo/solar_system/src/solar_system.h +++ b/demo/solar_system/src/solar_system.h @@ -25,6 +25,7 @@ namespace astrophysics { Real3 RotateVector(const Real3& vec, real_t th, char axis); inline int Simulate(int argc, const char** argv) { + constexpr real_t kSolarMassKg = 1.9884e30; Simulation simulation(argc, argv); auto* rm = simulation.GetResourceManager(); @@ -38,7 +39,7 @@ inline int Simulate(int argc, const char** argv) { // create sun Star* sun = new Star(109); sun->SetPosition({0, 0, 0}); - sun->SetMass(1988400e+24); + sun->SetMass(kSolarMassKg); sun->AddBehavior(new Gravity); rm->AddAgent(sun); diff --git a/demo/soma_clustering/src/my_cell.h b/demo/soma_clustering/src/my_cell.h index ada73ba66..b3fa1c89b 100644 --- a/demo/soma_clustering/src/my_cell.h +++ b/demo/soma_clustering/src/my_cell.h @@ -23,7 +23,7 @@ namespace soma_clustering { // Define my custom cell, which extends Cell by adding an extra // data member cell_type. class MyCell : public Cell { - BDM_AGENT_HEADER(MyCell, Cell, 1); + BDM_AGENT_HEADER(MyCell, Cell); public: MyCell() {} @@ -33,6 +33,15 @@ class MyCell : public Cell { int GetCellType() const { return cell_type_; } + bool GetVisualizationData(const std::string& name, + VisualizationData* values) const override { + if (name == "cell_type_") { + *values = std::vector{cell_type_}; + return true; + } + return Base::GetVisualizationData(name, values); + } + private: int cell_type_; }; diff --git a/demo/tumor_concept/src/tumor_concept.h b/demo/tumor_concept/src/tumor_concept.h index 7a3ec35be..012d9cb6c 100644 --- a/demo/tumor_concept/src/tumor_concept.h +++ b/demo/tumor_concept/src/tumor_concept.h @@ -27,7 +27,7 @@ namespace tumor_concept { // members: cell_color and can_divide class MyCell : public Cell { // our object extends the Cell object // create the header with our new data member - BDM_AGENT_HEADER(MyCell, Cell, 1); + BDM_AGENT_HEADER(MyCell, Cell); public: MyCell() {} @@ -56,6 +56,15 @@ class MyCell : public Cell { // our object extends the Cell object void SetCellColor(int cell_color) { cell_color_ = cell_color; } int GetCellColor() const { return cell_color_; } + bool GetVisualizationData(const std::string& name, + VisualizationData* values) const override { + if (name == "cell_color_") { + *values = std::vector{cell_color_}; + return true; + } + return Base::GetVisualizationData(name, values); + } + private: // declare new data member and define their type // private data can only be accessed by public function and not directly @@ -65,7 +74,7 @@ class MyCell : public Cell { // our object extends the Cell object // Define growth behaviour struct Growth : public Behavior { - BDM_BEHAVIOR_HEADER(Growth, Behavior, 1); + BDM_BEHAVIOR_HEADER(Growth, Behavior); Growth() { AlwaysCopyToNew(); } virtual ~Growth() {} diff --git a/doc/user_guide/tumor_concept.md b/doc/user_guide/tumor_concept.md index 80a5025f1..8f0282a3e 100644 --- a/doc/user_guide/tumor_concept.md +++ b/doc/user_guide/tumor_concept.md @@ -255,6 +255,15 @@ class MyCell : public Cell { // our object extends the Cell object void SetCellColor(int cell_color) { cell_color_ = cell_color; } int GetCellColor() const { return cell_color_; } + bool GetVisualizationData(const std::string& name, + VisualizationData* values) const override { + if (name == "cell_color_") { + *values = std::vector{cell_color_}; + return true; + } + return Base::GetVisualizationData(name, values); + } + private: // declare new data member and define their type // private data can only be accessed by public function and not directly diff --git a/src/core/agent/agent.cc b/src/core/agent/agent.cc index b2ec52bc5..5b0f401d9 100644 --- a/src/core/agent/agent.cc +++ b/src/core/agent/agent.cc @@ -33,7 +33,6 @@ #include "core/simulation.h" #include "core/util/log.h" #include "core/util/macros.h" -#include "core/util/root.h" #include "core/util/type.h" namespace bdm { @@ -42,8 +41,6 @@ Agent::Agent() { uid_ = Simulation::GetActive()->GetAgentUidGenerator()->GenerateUid(); } -Agent::Agent(TRootIOCtor* io_ctor) {} - Agent::Agent(const Agent& other) : uid_(other.uid_), box_idx_(other.box_idx_), @@ -122,6 +119,24 @@ void Agent::AssignNewUid() { const AgentUid& Agent::GetUid() const { return uid_; } +bool Agent::GetVisualizationData(const std::string& name, + VisualizationData* values) const { + if (name == "position_") { + const auto& position = GetPosition(); + *values = std::vector(position.begin(), position.end()); + return true; + } + if (name == "diameter_") { + *values = std::vector{GetDiameter()}; + return true; + } + if (name == "uid_") { + *values = std::vector{GetUid()}; + return true; + } + return false; +} + uint32_t Agent::GetBoxIdx() const { return box_idx_; } void Agent::SetBoxIdx(uint32_t idx) { box_idx_ = idx; } diff --git a/src/core/agent/agent.h b/src/core/agent/agent.h index 8b6975e84..067354641 100644 --- a/src/core/agent/agent.h +++ b/src/core/agent/agent.h @@ -24,6 +24,7 @@ #include #include #include +#include #include #include "core/agent/agent_pointer.h" @@ -34,7 +35,6 @@ #include "core/interaction_force.h" #include "core/shape.h" #include "core/util/macros.h" -#include "core/util/root.h" #include "core/util/spinlock.h" #include "core/util/type.h" @@ -43,18 +43,9 @@ namespace bdm { /// Macro to insert required boilerplate code into agent /// @param class_name scalar class name of the agent /// @param base_class scalar class name of the base agent -/// @param class_version_id required for ROOT I/O (see ROOT BDM_CLASS_DEF -/// Macro). -/// Every time the layout of the class is changed, class_version_id -/// must be incremented by one. The class_version_id should be greater -/// or equal to 1. -/// @param ...: List of all data members of this class -#define BDM_AGENT_HEADER(class_name, base_class, class_version_id) \ +#define BDM_AGENT_HEADER(class_name, base_class) \ public: \ using Base = base_class; \ - \ - explicit class_name(TRootIOCtor* io_ctor) {} \ - \ /** Create a new instance of this object using the default constructor. */ \ Agent* New() const override { return new class_name(); } \ /** Create a new instance of this object using the copy constructor. */ \ @@ -67,21 +58,20 @@ namespace bdm { Base* UpCast() { return static_cast(this); } \ \ /** Cast `this` to the base class pointer (one level up) */ \ - const Base* UpCast() const { return static_cast(this); } \ - \ - BDM_CLASS_DEF_OVERRIDE(class_name, class_version_id) + const Base* UpCast() const { return static_cast(this); } // ----------------------------------------------------------------------------- class Behavior; +using VisualizationData = + std::variant, std::vector, std::vector>; + /// Contains code required by all agents class Agent { public: Agent(); - explicit Agent(TRootIOCtor* io_ctor); - Agent(const Agent& other); virtual ~Agent(); @@ -141,6 +131,10 @@ class Agent { return {"position_", "diameter_"}; } + /// Copies one typed visualization value into `values`. + virtual bool GetVisualizationData(const std::string& name, + VisualizationData* values) const; + virtual void RunDiscretization(); void AssignNewUid(); @@ -287,19 +281,19 @@ class Agent { } private: - Spinlock lock_; //! + Spinlock lock_; /// Helper variable used to support removal of behaviors while /// `RunBehaviors` iterates over them. uint16_t run_behavior_loop_idx_ = 0; /// If an agent is static, we should not compute the mechanical forces - bool is_static_ = false; //! + bool is_static_ = false; /// If an agent becomes non-static (i.e. it moved or grew), we should set this /// flag to true to also compute mechanical forces on the neighboring agents - bool propagate_staticness_neighborhood_ = true; //! + bool propagate_staticness_neighborhood_ = true; /// Flag to determine of an agent is static in the next timestep - mutable bool is_static_next_ts_ = false; //! + mutable bool is_static_next_ts_ = false; /// Function to copy behaviors from existing Agent to this one /// and to initialize them. @@ -314,8 +308,6 @@ class Agent { /// This function sets the attributes `NewAgentEvent::existing_behavior` /// and `NewAgentEvent::new_behaviors` to their correct value. void UpdateBehaviors(const NewAgentEvent& event); - - BDM_CLASS_DEF(Agent, 1) }; } // namespace bdm diff --git a/src/core/agent/agent_handle.h b/src/core/agent/agent_handle.h index d15971ea4..fd8fac454 100644 --- a/src/core/agent/agent_handle.h +++ b/src/core/agent/agent_handle.h @@ -15,8 +15,9 @@ #ifndef CORE_AGENT_AGENT_HANDLE_H_ #define CORE_AGENT_AGENT_HANDLE_H_ +#include #include -#include "core/util/root.h" +#include namespace bdm { @@ -72,8 +73,6 @@ class AgentHandle { /// changed element index to uint32_t after issues with std::atomic with /// size 16 -> max element_idx: 4.294.967.296 ElementIdx_t element_idx_; - - BDM_CLASS_DEF_NV(AgentHandle, 1); }; } // namespace bdm diff --git a/src/core/agent/agent_pointer.h b/src/core/agent/agent_pointer.h index f94a198a3..ee39eae5a 100644 --- a/src/core/agent/agent_pointer.h +++ b/src/core/agent/agent_pointer.h @@ -23,7 +23,6 @@ #include "core/agent/agent_uid.h" #include "core/execution_context/execution_context.h" #include "core/simulation.h" -#include "core/util/root.h" namespace bdm { @@ -240,7 +239,7 @@ class AgentPointer { }; private: - Data d_ = {AgentUid()}; //! + Data d_ = {AgentUid()}; template typename std::enable_if::value, TTo*>::type Cast( @@ -253,43 +252,8 @@ class AgentPointer { TFrom* agent) const { return dynamic_cast(agent); } - - BDM_CLASS_DEF_NV(AgentPointer, 3); }; -// The following custom streamer should be visible to rootcling for dictionary -// generation, but not to the interpreter! -#if (!defined(__CLING__) || defined(__ROOTCLING__)) && defined(USE_DICT) - -template -inline void AgentPointer::Streamer(TBuffer& R__b) { - if (R__b.IsReading()) { - R__b.ReadClassBuffer(AgentPointer::Class(), this); - AgentUid restored_uid; - R__b.ReadClassBuffer(AgentUid::Class(), &restored_uid); - if (gAgentPointerMode == AgentPointerMode::kIndirect) { - d_.uid = restored_uid; - } else if (restored_uid != AgentUid()) { - auto* ctxt = Simulation::GetActive()->GetExecutionContext(); - d_.agent = Cast(ctxt->GetAgent(restored_uid)); - } else { - d_.agent = nullptr; - } - } else { - R__b.WriteClassBuffer(AgentPointer::Class(), this); - AgentUid uid; - if (gAgentPointerMode == AgentPointerMode::kIndirect) { - uid = d_.uid; - R__b.WriteClassBuffer(AgentUid::Class(), &d_.uid); - } else if (d_.agent != nullptr) { - uid = d_.agent->GetUid(); - } - R__b.WriteClassBuffer(AgentUid::Class(), &uid); - } -} - -#endif // !defined(__CLING__) || defined(__ROOTCLING__) - template struct is_agent_ptr { // NOLINT static constexpr bool value = false; // NOLINT diff --git a/src/core/agent/agent_uid.h b/src/core/agent/agent_uid.h index 67f86946d..b89d3edd5 100644 --- a/src/core/agent/agent_uid.h +++ b/src/core/agent/agent_uid.h @@ -15,8 +15,10 @@ #ifndef CORE_AGENT_AGENT_UID_H_ #define CORE_AGENT_AGENT_UID_H_ +#include +#include #include -#include "core/util/root.h" +#include namespace bdm { @@ -108,8 +110,6 @@ class AgentUid { /// Determines how often index_ has been resused Reused_t reused_; - - BDM_CLASS_DEF_NV(AgentUid, 1); }; } // namespace bdm diff --git a/src/core/agent/agent_uid_generator.h b/src/core/agent/agent_uid_generator.h index 96c2dfd20..38a5f6f22 100644 --- a/src/core/agent/agent_uid_generator.h +++ b/src/core/agent/agent_uid_generator.h @@ -24,7 +24,6 @@ #include "core/container/shared_data.h" #include "core/scheduler.h" #include "core/simulation.h" -#include "core/util/root.h" #include "core/util/spinlock.h" namespace bdm { @@ -67,34 +66,13 @@ class AgentUidGenerator { void Update() { tl_uids_.resize(tinfo_->GetMaxThreads()); } private: - std::atomic counter_; //! - /// ROOT can't persist std::atomic. - /// Therefore this additional helper variable is needed. - typename AgentUid::Index_t root_counter_; + std::atomic counter_; /// Thread local vector of AgentUids that can be reused SharedData> tl_uids_; - ThreadInfo* tinfo_ = nullptr; //! - - BDM_CLASS_DEF_NV(AgentUidGenerator, 1); + ThreadInfo* tinfo_ = nullptr; }; -// The following custom streamer should be visible to rootcling for dictionary -// generation, but not to the interpreter! -#if (!defined(__CLING__) || defined(__ROOTCLING__)) && defined(USE_DICT) - -inline void AgentUidGenerator::Streamer(TBuffer& R__b) { - if (R__b.IsReading()) { - R__b.ReadClassBuffer(AgentUidGenerator::Class(), this); - this->counter_ = this->root_counter_; - } else { - this->root_counter_ = this->counter_.load(); - R__b.WriteClassBuffer(AgentUidGenerator::Class(), this); - } -} - -#endif // !defined(__CLING__) || defined(__ROOTCLING__) - } // namespace bdm #endif // CORE_AGENT_AGENT_UID_GENERATOR_H_ diff --git a/src/core/agent/cell.h b/src/core/agent/cell.h index 5f7407dd4..723d9889f 100644 --- a/src/core/agent/cell.h +++ b/src/core/agent/cell.h @@ -38,7 +38,7 @@ namespace bdm { class Cell : public Agent { - BDM_AGENT_HEADER(Cell, Agent, 1); + BDM_AGENT_HEADER(Cell, Agent); public: /// First axis of the local coordinate system. @@ -191,6 +191,28 @@ class Cell : public Agent { real_t GetVolume() const { return volume_; } + bool GetVisualizationData(const std::string& name, + VisualizationData* values) const override { + if (name == "volume_") { + *values = std::vector{volume_}; + return true; + } + if (name == "density_") { + *values = std::vector{density_}; + return true; + } + if (name == "adherence_") { + *values = std::vector{adherence_}; + return true; + } + if (name == "tractor_force_") { + *values = + std::vector(tractor_force_.begin(), tractor_force_.end()); + return true; + } + return Base::GetVisualizationData(name, values); + } + void SetAdherence(real_t adherence) { if (adherence < adherence_) { SetStaticnessNextTimestep(false); diff --git a/src/core/agent/spherical_agent.h b/src/core/agent/spherical_agent.h index 3bd125cfc..5491c3e4c 100644 --- a/src/core/agent/spherical_agent.h +++ b/src/core/agent/spherical_agent.h @@ -24,7 +24,7 @@ namespace bdm { class SphericalAgent : public Agent { - BDM_AGENT_HEADER(SphericalAgent, Agent, 1); + BDM_AGENT_HEADER(SphericalAgent, Agent); public: SphericalAgent() : diameter_(1.0) {} diff --git a/src/core/analysis/reduce.h b/src/core/analysis/reduce.h index 37ce6309c..4e776f532 100644 --- a/src/core/analysis/reduce.h +++ b/src/core/analysis/reduce.h @@ -41,7 +41,6 @@ struct Reducer : public Functor { /// Resets the internal state between calculations. virtual void Reset() = 0; virtual Reducer* NewCopy() const = 0; - BDM_CLASS_DEF(Reducer, 1) }; // ----------------------------------------------------------------------------- @@ -121,49 +120,13 @@ class GenericReducer : public Reducer { } private: - SharedData tl_results_; //! - void (*agent_function_)(Agent*, T*) = nullptr; //! - T (*reduce_partial_results_)(const SharedData&) = nullptr; //! - bool (*filter_)(Agent*) = nullptr; //! - TResult (*post_process_)(TResult) = nullptr; //! - BDM_CLASS_DEF_OVERRIDE(GenericReducer, 1) + SharedData tl_results_; + void (*agent_function_)(Agent*, T*) = nullptr; + T (*reduce_partial_results_)(const SharedData&) = nullptr; + bool (*filter_)(Agent*) = nullptr; + TResult (*post_process_)(TResult) = nullptr; }; -// The following custom streamer should be visible to rootcling for dictionary -// generation, but not to the interpreter! -#if (!defined(__CLING__) || defined(__ROOTCLING__)) && defined(USE_DICT) - -// The custom streamer is needed because ROOT can't stream function pointers -// by default. -template -inline void GenericReducer::Streamer(TBuffer& R__b) { - if (R__b.IsReading()) { - R__b.ReadClassBuffer(GenericReducer::Class(), this); - Long64_t l; - R__b.ReadLong64(l); - this->agent_function_ = reinterpret_cast(l); - R__b.ReadLong64(l); - this->reduce_partial_results_ = - reinterpret_cast&)>(l); - R__b.ReadLong64(l); - this->filter_ = reinterpret_cast(l); - R__b.ReadLong64(l); - this->post_process_ = reinterpret_cast(l); - } else { - R__b.WriteClassBuffer(GenericReducer::Class(), this); - Long64_t l = reinterpret_cast(this->agent_function_); - R__b.WriteLong64(l); - l = reinterpret_cast(this->reduce_partial_results_); - R__b.WriteLong64(l); - l = reinterpret_cast(this->filter_); - R__b.WriteLong64(l); - l = reinterpret_cast(this->post_process_); - R__b.WriteLong64(l); - } -} - -#endif // !defined(__CLING__) || defined(__ROOTCLING__) - /// Iterates over all agents executing the `agent_functor` and updating a /// a thread-local and therefore partial result. /// The second parameter specifies how these partial results should be combined @@ -270,38 +233,11 @@ struct Counter : public Reducer { Reducer* NewCopy() const override { return new Counter(*this); } private: - SharedData tl_results_; //! - bool (*condition_)(Agent*) = nullptr; //! - TResult (*post_process_)(TResult) = nullptr; //! - BDM_CLASS_DEF_OVERRIDE(Counter, 1) + SharedData tl_results_; + bool (*condition_)(Agent*) = nullptr; + TResult (*post_process_)(TResult) = nullptr; }; -// The following custom streamer should be visible to rootcling for dictionary -// generation, but not to the interpreter! -#if (!defined(__CLING__) || defined(__ROOTCLING__)) && defined(USE_DICT) - -// The custom streamer is needed because ROOT can't stream function pointers -// by default. -template -inline void Counter::Streamer(TBuffer& R__b) { - if (R__b.IsReading()) { - R__b.ReadClassBuffer(Counter::Class(), this); - Long64_t l; - R__b.ReadLong64(l); - this->condition_ = reinterpret_cast(l); - R__b.ReadLong64(l); - this->post_process_ = reinterpret_cast(l); - } else { - R__b.WriteClassBuffer(Counter::Class(), this); - Long64_t l = reinterpret_cast(this->condition_); - R__b.WriteLong64(l); - l = reinterpret_cast(this->post_process_); - R__b.WriteLong64(l); - } -} - -#endif // !defined(__CLING__) || defined(__ROOTCLING__) - /// Counts the number of agents for which `condition` evaluates to true. /// Let's assume we want to count all infected agents in a virus spreading /// simulation. diff --git a/src/core/behavior/behavior.h b/src/core/behavior/behavior.h index 8a9dcc738..fd2bd1e12 100644 --- a/src/core/behavior/behavior.h +++ b/src/core/behavior/behavior.h @@ -123,22 +123,16 @@ class Behavior { private: NewAgentEventUid copy_mask_ = 0; NewAgentEventUid remove_mask_ = 0; - BDM_CLASS_DEF(Behavior, 2); }; /// Inserts boilerplate code for behaviors with state -#define BDM_BEHAVIOR_HEADER(class_name, base_class, class_version_id) \ +#define BDM_BEHAVIOR_HEADER(class_name, base_class) \ public: \ using Base = base_class; \ /** Create a new instance of this object using the default constructor. */ \ Behavior* New() const override { return new class_name(); } \ /** Create a new instance of this object using the copy constructor. */ \ - Behavior* NewCopy() const override { return new class_name(*this); } \ - \ - private: \ - BDM_CLASS_DEF_OVERRIDE(class_name, class_version_id); \ - \ - public: + Behavior* NewCopy() const override { return new class_name(*this); } } // namespace bdm diff --git a/src/core/behavior/chemotaxis.h b/src/core/behavior/chemotaxis.h index ca966b273..864af07a2 100644 --- a/src/core/behavior/chemotaxis.h +++ b/src/core/behavior/chemotaxis.h @@ -24,7 +24,7 @@ namespace bdm { /// Move cells along the diffusion gradient (from low concentration to high) class Chemotaxis : public Behavior { - BDM_BEHAVIOR_HEADER(Chemotaxis, Behavior, 1); + BDM_BEHAVIOR_HEADER(Chemotaxis, Behavior); public: Chemotaxis() = default; diff --git a/src/core/behavior/gene_regulation.h b/src/core/behavior/gene_regulation.h index 59879257a..ed114f0b7 100644 --- a/src/core/behavior/gene_regulation.h +++ b/src/core/behavior/gene_regulation.h @@ -22,7 +22,6 @@ #include "core/param/param.h" #include "core/scheduler.h" #include "core/simulation.h" -#include "core/util/root.h" namespace bdm { @@ -34,7 +33,7 @@ namespace bdm { /// The user determines which method is picked in particular simulation /// through variable `Param::numerical_ode_solver`. class GeneRegulation : public Behavior { - BDM_BEHAVIOR_HEADER(GeneRegulation, Behavior, 1); + BDM_BEHAVIOR_HEADER(GeneRegulation, Behavior); public: GeneRegulation() { AlwaysCopyToNew(); } diff --git a/src/core/behavior/growth_division.h b/src/core/behavior/growth_division.h index d82acdd14..d23a3fe1e 100644 --- a/src/core/behavior/growth_division.h +++ b/src/core/behavior/growth_division.h @@ -19,14 +19,13 @@ #include "core/agent/cell_division_event.h" #include "core/behavior/behavior.h" #include "core/util/log.h" -#include "core/util/root.h" namespace bdm { /// This behavior grows the agent until the diameter reaches /// the specified threshold and divides the object afterwards. class GrowthDivision : public Behavior { - BDM_BEHAVIOR_HEADER(GrowthDivision, Behavior, 1); + BDM_BEHAVIOR_HEADER(GrowthDivision, Behavior); public: GrowthDivision() { AlwaysCopyToNew(); } diff --git a/src/core/behavior/regulatory_network.h b/src/core/behavior/regulatory_network.h index 1914ce29e..fe0a6b37c 100644 --- a/src/core/behavior/regulatory_network.h +++ b/src/core/behavior/regulatory_network.h @@ -19,27 +19,22 @@ #include "core/behavior/behavior.h" -#ifndef __ROOTCLING__ #include "boost/numeric/odeint.hpp" #include "boost/phoenix/core.hpp" #include "boost/phoenix/operator.hpp" -#endif -#ifndef __ROOTCLING__ typedef boost::numeric::ublas::vector b_vector_t; typedef boost::numeric::ublas::matrix b_matrix_t; -#endif namespace bdm { enum class ODE_solver { Euler, Rosenbrock, RungeKutta }; class RegulatoryNetwork : public Behavior { - BDM_BEHAVIOR_HEADER(RegulatoryNetwork, Behavior, 1); + BDM_BEHAVIOR_HEADER(RegulatoryNetwork, Behavior); public: RegulatoryNetwork() { AlwaysCopyToNew(); } -#ifndef __ROOTCLING__ RegulatoryNetwork( real_t dt, int n_dt, const std::vector& x, ODE_solver m, const std::function& @@ -58,7 +53,6 @@ class RegulatoryNetwork : public Behavior { out_ = out; method_ = m; } -#endif virtual ~RegulatoryNetwork() = default; void Initialize(const NewAgentEvent& event) override { @@ -66,35 +60,28 @@ class RegulatoryNetwork : public Behavior { if (auto* other = dynamic_cast(event.existing_behavior)) { -#ifndef __ROOTCLING__ current_time_ = other->current_time_; current_species_ = other->current_species_; previous_species_ = other->previous_species_; -#endif time_step_ = other->time_step_; time_subdivision_ = other->time_subdivision_; -#ifndef __ROOTCLING__ rhs_ = other->rhs_; jacob_ = other->jacob_; out_ = other->out_; method_ = other->method_; -#endif } else { Log::Fatal("RegulatoryNetwork::EventConstructor", "other was not of type RegulatoryNetwork"); } } -#ifndef __ROOTCLING__ const size_t GetNumberOfSpecies() const { return current_species_.size(); } const b_vector_t& GetSpecies() const { return current_species_; } const real_t& GetSpecie(size_t i) const { return current_species_[i]; } -#endif void Run(Agent* agent) override { -#ifndef __ROOTCLING__ // update the previous solution previous_species_ = current_species_; @@ -153,14 +140,9 @@ class RegulatoryNetwork : public Behavior { // print-out the results out_(current_species_, current_time_, agent); -#else - Log::Fatal("RegulatoryNetwork::Run", - "this behavior is supported only with \"boost\" installed"); -#endif }; protected: -#ifndef __ROOTCLING__ void SetInitialSpecies(const std::vector& x) { const size_t n_species = x.size(); @@ -169,7 +151,6 @@ class RegulatoryNetwork : public Behavior { for (size_t i = 0; i < n_species; i++) current_species_[i] = previous_species_[i] = x[i]; } -#endif private: /// Pseudo-time for ODE(s) time integration @@ -177,22 +158,18 @@ class RegulatoryNetwork : public Behavior { /// Time-step for ODE(s) time integration real_t time_step_ = 1.0; int time_subdivision_ = 100; -#ifndef __ROOTCLING__ /// Current solution of the species concentration b_vector_t current_species_ = {}; /// Previous solution of the species concentration b_vector_t previous_species_ = {}; /// Method used for the ODE(s) numerical solution ODE_solver method_; -#endif -#ifndef __ROOTCLING__ std::function rhs_; std::function jacob_; std::function out_; -#endif }; } // namespace bdm diff --git a/src/core/behavior/secretion.h b/src/core/behavior/secretion.h index 50813c1b3..c7e89ee90 100644 --- a/src/core/behavior/secretion.h +++ b/src/core/behavior/secretion.h @@ -26,7 +26,7 @@ namespace bdm { /// Secrete substance at Agent position class Secretion : public Behavior { - BDM_BEHAVIOR_HEADER(Secretion, Behavior, 2); + BDM_BEHAVIOR_HEADER(Secretion, Behavior); public: Secretion() = default; diff --git a/src/core/behavior/stateless_behavior.h b/src/core/behavior/stateless_behavior.h index 4fb32ad35..3bcc3fe05 100644 --- a/src/core/behavior/stateless_behavior.h +++ b/src/core/behavior/stateless_behavior.h @@ -33,7 +33,7 @@ namespace bdm { /// Without StatelessBehavior the following code would be required. /// \code /// struct RapidDivision : public Behavior { -/// BDM_BEHAVIOR_HEADER(RapidDivision, Behavior, 1); +/// BDM_BEHAVIOR_HEADER(RapidDivision, Behavior); /// /// RapidDivision() = default; /// virtual ~RapidDivision() = default; @@ -44,7 +44,7 @@ namespace bdm { /// }; /// \endcode class StatelessBehavior : public Behavior { - BDM_BEHAVIOR_HEADER(StatelessBehavior, Behavior, 1); + BDM_BEHAVIOR_HEADER(StatelessBehavior, Behavior); public: using FPtr = void (*)(Agent*); @@ -67,30 +67,9 @@ class StatelessBehavior : public Behavior { } private: - FPtr fptr_; //! + FPtr fptr_; }; -// The following custom streamer should be visible to rootcling for dictionary -// generation, but not to the interpreter! -#if (!defined(__CLING__) || defined(__ROOTCLING__)) && defined(USE_DICT) - -// The custom streamer is needed because ROOT can't stream function pointers -// by default. -inline void StatelessBehavior::Streamer(TBuffer& R__b) { - if (R__b.IsReading()) { - R__b.ReadClassBuffer(StatelessBehavior::Class(), this); - Long64_t l; - R__b.ReadLong64(l); - this->fptr_ = reinterpret_cast(l); - } else { - R__b.WriteClassBuffer(StatelessBehavior::Class(), this); - Long64_t l = reinterpret_cast(this->fptr_); - R__b.WriteLong64(l); - } -} - -#endif // !defined(__CLING__) || defined(__ROOTCLING__) - } // namespace bdm #endif // CORE_BEHAVIOR_STATELESS_BEHAVIOR_H_ diff --git a/src/core/model_initializer.h b/src/core/model_initializer.h index 578c2b4c9..53947fcc8 100644 --- a/src/core/model_initializer.h +++ b/src/core/model_initializer.h @@ -15,12 +15,12 @@ #ifndef CORE_MODEL_INITIALIZER_H_ #define CORE_MODEL_INITIALIZER_H_ -#include #include #include #include #include +#include "core/container/fixed_size_vector.h" #include "core/container/math_array.h" #include "core/diffusion/diffusion_grid.h" #include "core/resource_manager.h" @@ -335,42 +335,14 @@ struct ModelInitializer { static void CreateAgentsInSphereRndm(const Real3& center, real_t radius, uint64_t num_agents, Function agent_builder) { - // We use a probability density function (PDF) to model the probability of - // an agent to occur at a distance `r>=0` of the center. As the surface of - // a sphere scales as `r^2`, the PDF does as well. Thus - // `p(r)=a*r^2*\Theta(R-r)`, where `\Theta` is a heavyside function and R is - // largest allowed radius (interpretation: no agents outside the sphere). We - // can fix `a` by requiring `\int_0^\inf p(r') dr' = 1` and obtain - // `a=3/R^3`. - auto radial_pdf_sphere = [](const double* x, const double* params) { - double R{params[0]}; - double r{x[0]}; - if (r > 0.0 && r <= R) { - return 3.0 * std::pow(r, 2.0) / std::pow(R, 3.0); - } else { - return 0.0; - } - }; - - // Get a random number generator to sample from our PDF. - auto* random = Simulation::GetActive()->GetRandom(); - auto rng = - random->GetUserDefinedDistRng1D(radial_pdf_sphere, {radius}, 0, radius); - - // Create a random radius for each of the agents. Note: this is done - // serially because we GetUserDefinedDistRng1D does not work in parallel - // regions at the moment. - std::vector random_radius; - random_radius.resize(num_agents); - for (size_t i = 0; i < num_agents; i++) { - random_radius[i] = rng.Sample(); - } -#pragma omp parallel shared(random_radius) +#pragma omp parallel { + auto* random = Simulation::GetActive()->GetRandom(); auto* ctxt_tl = Simulation::GetActive()->GetExecutionContext(); #pragma omp for schedule(static) for (uint64_t i = 0; i < num_agents; i++) { - auto pos = random->Sphere(random_radius[i]) + center; + auto radial_distance = radius * std::cbrt(random->Uniform()); + auto pos = random->Sphere(radial_distance) + center; auto* new_agent = agent_builder(pos); ctxt_tl->AddAgent(new_agent); } diff --git a/src/core/type_index.cc b/src/core/type_index.cc index 7ff1b4db9..7077d9848 100644 --- a/src/core/type_index.cc +++ b/src/core/type_index.cc @@ -14,13 +14,11 @@ #include "core/type_index.h" -#include - namespace bdm { // ----------------------------------------------------------------------------- void TypeIndex::Add(Agent* agent) { - auto& type_vector = data_[agent->IsA()]; + auto& type_vector = data_[agent->GetTypeName()]; auto uid = agent->GetUid(); if (index_.size() <= uid.GetIndex()) { Reserve(uid.GetIndex() + 1); @@ -32,14 +30,14 @@ void TypeIndex::Add(Agent* agent) { // ----------------------------------------------------------------------------- void TypeIndex::Update(Agent* new_agent) { auto idx = index_[new_agent->GetUid()]; - auto& type_vector = data_[new_agent->IsA()]; + auto& type_vector = data_[new_agent->GetTypeName()]; type_vector[idx] = new_agent; } // ----------------------------------------------------------------------------- void TypeIndex::Remove(Agent* agent) { auto idx = index_[agent->GetUid()]; - auto& type_vector = data_[agent->IsA()]; + auto& type_vector = data_[agent->GetTypeName()]; if (idx == type_vector.size() - 1) { type_vector.pop_back(); } else { @@ -67,8 +65,9 @@ void TypeIndex::Reserve(uint64_t capacity) { } // ----------------------------------------------------------------------------- -const std::vector& TypeIndex::GetType(TClass* tclass) const { - return data_[tclass]; +const std::vector& TypeIndex::GetType( + const std::string& type_name) const { + return data_[type_name]; } } // namespace bdm diff --git a/src/core/type_index.h b/src/core/type_index.h index 5c611a1ce..c90f68181 100644 --- a/src/core/type_index.h +++ b/src/core/type_index.h @@ -15,13 +15,12 @@ #ifndef CORE_TYPE_INDEX_H_ #define CORE_TYPE_INDEX_H_ +#include #include #include "core/agent/agent.h" #include "core/container/agent_uid_map.h" #include "core/container/flatmap.h" -class TClass; - namespace bdm { class TypeIndex { @@ -36,10 +35,10 @@ class TypeIndex { void Reserve(uint64_t capacity); - const std::vector& GetType(TClass* tclass) const; + const std::vector& GetType(const std::string& type_name) const; private: - UnorderedFlatmap> data_; + UnorderedFlatmap> data_; AgentUidMap index_; }; diff --git a/src/core/util/log.h b/src/core/util/log.h index 4d4a17a6e..e2c699514 100644 --- a/src/core/util/log.h +++ b/src/core/util/log.h @@ -5,7 +5,6 @@ // // Licensed under the Apache License, Version 2.0 (the "License"); // you may not use this file except in compliance with the License. -// // See the LICENSE file distributed with this work for details. // See the NOTICE file distributed with this work for additional information // regarding copyright ownership. @@ -15,122 +14,86 @@ #ifndef CORE_UTIL_LOG_H_ #define CORE_UTIL_LOG_H_ -#include - +#include +#include #include #include +#include #include -#include #include #include "core/util/string.h" namespace bdm { -/// @brief Wrapper class over ROOT logging module -/// class Log { public: - /// @brief Prints debug message - /// - /// @param[in] location The location of the message - /// @param[in] parts objects that compose the entire message - /// + enum class Level { kDebug, kInfo, kWarning, kError }; + + static void SetLevel(Level level) { level_.store(level); } + static Level GetLevel() { return level_.load(); } + template static void Debug(const std::string& location, const Args&... parts) { - // kPrint has the highest level of verbosity - if (gErrorIgnoreLevel <= kPrint) { - std::string message = Concat(parts...); - // Mimic ROOT logging output - fprintf(stderr, "Debug in <%s>: %s\n", location.c_str(), message.c_str()); - } + Write(Level::kDebug, "Debug", location, Concat(parts...)); } - /// @brief Prints information message - /// - /// @param[in] location The location of the message - /// @param[in] parts objects that compose the entire message - /// template static void Info(const std::string& location, const Args&... parts) { - std::string message = Concat(parts...); - // ROOT function - ::Info(location.c_str(), "%s", message.c_str()); + Write(Level::kInfo, "Info", location, Concat(parts...)); } - /// @brief Prints warning message - /// - /// @param[in] location The location of the message - /// @param[in] parts objects that compose the entire message - /// template static void Warning(const std::string& location, const Args&... parts) { - std::string message = Concat(parts...); - // ROOT function - ::Warning(location.c_str(), "%s", message.c_str()); + Write(Level::kWarning, "Warning", location, Concat(parts...)); } - /// @brief Prints error message - /// - /// @param[in] location The location of the message - /// @param[in] parts objects that compose the entire message - /// template static void Error(const std::string& location, const Args&... parts) { - std::string message = Concat(parts...); - // ROOT function - ::Error(location.c_str(), "%s", message.c_str()); + Write(Level::kError, "Error", location, Concat(parts...)); } - /// @brief Prints break message - /// - /// @param[in] location The location of the message - /// @param[in] parts objects that compose the entire message - /// template static void Break(const std::string& location, const Args&... parts) { - std::string message = Concat(parts...); - // ROOT function - ::Break(location.c_str(), "%s", message.c_str()); + Write(Level::kError, "Break", location, Concat(parts...)); } - /// @brief Prints system error message - /// - /// @param[in] location The location of the message - /// @param[in] parts objects that compose the entire message - /// template static void SysError(const std::string& location, const Args&... parts) { - std::string message = Concat(parts...); - // ROOT function - ::SysError(location.c_str(), "%s", message.c_str()); + Write(Level::kError, "SysError", location, + Concat(parts..., ": ", std::strerror(errno))); } - /// @brief Prints fatal error message - /// - /// @param[in] location The location of the message - /// @param[in] parts objects that compose the entire message - /// template - static void Fatal(const std::string& location, const Args&... parts) { + [[noreturn]] static void Fatal(const std::string& location, + const Args&... parts) { std::string message = Concat(parts...); - // ROOT function - ::Error(location.c_str(), "%s", message.c_str()); - // std::runtime_error() will fail the DeathTests - exit(1); + std::fprintf(stderr, "Error in <%s>: %s\n", location.c_str(), + message.c_str()); + std::exit(1); } template static void Condition(const std::function& lambda, const Args&... parts) { - // kPrint has the highest level of verbosity if (lambda()) { std::string message = Concat(parts...); - // Mimic ROOT logging output - fprintf(stdout, "%s\n", message.c_str()); + std::fprintf(stdout, "%s\n", message.c_str()); + } + } + + private: + static void Write(Level level, const char* name, const std::string& location, + const std::string& message) { + if (level >= level_.load()) { + std::fprintf(stderr, "%s in <%s>: %s\n", name, location.c_str(), + message.c_str()); } } + + inline static std::atomic level_ = Level::kWarning; }; + } // namespace bdm #endif // CORE_UTIL_LOG_H_ diff --git a/src/neuroscience/neurite_element.cc b/src/neuroscience/neurite_element.cc index 9ab8cc1b2..14a9443e2 100644 --- a/src/neuroscience/neurite_element.cc +++ b/src/neuroscience/neurite_element.cc @@ -15,6 +15,8 @@ #include "neuroscience/neurite_element.h" #include +#include "core/util/log.h" + namespace bdm { namespace neuroscience { @@ -135,6 +137,27 @@ std::set NeuriteElement::GetRequiredVisDataMembers() const { return {"mass_location_", "diameter_", "actual_length_", "spring_axis_"}; } +bool NeuriteElement::GetVisualizationData(const std::string& name, + VisualizationData* values) const { + if (name == "mass_location_") { + *values = std::vector(mass_location_.begin(), mass_location_.end()); + return true; + } + if (name == "actual_length_") { + *values = std::vector{actual_length_}; + return true; + } + if (name == "spring_axis_") { + *values = std::vector(spring_axis_.begin(), spring_axis_.end()); + return true; + } + if (name == "daughter_right_") { + *values = std::vector{daughter_right_.GetUidAsUint64()}; + return true; + } + return Base::GetVisualizationData(name, values); +} + void NeuriteElement::SetDiameter(real_t diameter) { if (diameter > diameter_) { SetPropagateStaticness(); @@ -285,8 +308,8 @@ std::array NeuriteElement::Bifurcate( // 1) physical bifurcation // check it is a terminal branch if (daughter_left_ != nullptr) { - Fatal("NeuriteElements", - "Bifurcation only allowed on a terminal neurite element"); + Log::Fatal("NeuriteElements", + "Bifurcation only allowed on a terminal neurite element"); } NeuriteBifurcationEvent event(length, diameter_1, diameter_2, direction_1, direction_2); @@ -347,7 +370,7 @@ void NeuriteElement::RemoveDaughter( daughter_right_ = nullptr; return; } - Fatal("NeuriteElement", "Given object is not a daughter!"); + Log::Fatal("NeuriteElement", "Given object is not a daughter!"); } void NeuriteElement::UpdateRelative(const NeuronOrNeurite& old_relative, @@ -369,7 +392,7 @@ void NeuriteElement::UpdateRelative(const NeuronOrNeurite& old_relative, Real3 NeuriteElement::ForceTransmittedFromDaugtherToMother( const NeuronOrNeurite& mother) { if (mother_ != &mother) { - Fatal("NeuriteElement", "Given object is not the mother!"); + Log::Fatal("NeuriteElement", "Given object is not the mother!"); return {0, 0, 0}; } @@ -903,8 +926,9 @@ void NeuriteElement::RemoveProximalNeuriteElement() { NeuriteElement* NeuriteElement::ExtendSideNeuriteElement( real_t length, real_t diameter, const Real3& direction) { if (daughter_right_ != nullptr) { - Fatal("NeuriteElement", - "Can't extend a side neurite since daughter_right is not a nullptr!"); + Log::Fatal( + "NeuriteElement", + "Can't extend a side neurite since daughter_right is not a nullptr!"); } SideNeuriteExtensionEvent event{length, diameter, direction}; diff --git a/src/neuroscience/neurite_element.h b/src/neuroscience/neurite_element.h index aafd9817a..bce49d165 100644 --- a/src/neuroscience/neurite_element.h +++ b/src/neuroscience/neurite_element.h @@ -54,7 +54,7 @@ namespace neuroscience { /// Only the distal end is moved. All the forces that are applied to the /// proximal node are transmitted to the mother element class NeuriteElement : public Agent, public NeuronOrNeurite { - BDM_AGENT_HEADER(NeuriteElement, Agent, 1); + BDM_AGENT_HEADER(NeuriteElement, Agent); public: NeuriteElement(); @@ -75,6 +75,9 @@ class NeuriteElement : public Agent, public NeuronOrNeurite { /// object. std::set GetRequiredVisDataMembers() const override; + bool GetVisualizationData(const std::string& name, + VisualizationData* values) const override; + void SetDiameter(real_t diameter) override; void SetDensity(real_t density); @@ -456,8 +459,6 @@ class NeuriteElement : public Agent, public NeuronOrNeurite { void Copy(const NeuriteElement& rhs); private: - // TODO(lukas) data members same as in cell -> resolve once ROOT-9321 has been - // resolved /// mass_location_ is distal end of the cylinder /// NB: Use setter and don't assign values directly Real3 mass_location_ = {{0.0, 0.0, 0.0}}; diff --git a/src/neuroscience/neuron_soma.cc b/src/neuroscience/neuron_soma.cc index b852e7012..58b1118c9 100644 --- a/src/neuroscience/neuron_soma.cc +++ b/src/neuroscience/neuron_soma.cc @@ -19,6 +19,7 @@ #include #include "core/resource_manager.h" +#include "core/util/log.h" #include "neuroscience/neurite_element.h" #include "neuroscience/new_agent_event/new_neurite_extension_event.h" #include "neuroscience/param.h" @@ -38,10 +39,11 @@ void NeuronSoma::Initialize(const NewAgentEvent& event) { if (event.GetUid() == CellDivisionEvent::kUid) { auto* mother = bdm_static_cast(event.existing_agent); if (mother->daughters_.size() != 0) { - Fatal("NeuronSoma", - "Dividing a neuron soma with attached neurites is not supported " - "in the default implementation! If you want to change this " - "behavior derive from this class and overwrite this method."); + Log::Fatal( + "NeuronSoma", + "Dividing a neuron soma with attached neurites is not supported " + "in the default implementation! If you want to change this " + "behavior derive from this class and overwrite this method."); } } } diff --git a/src/neuroscience/neuron_soma.h b/src/neuroscience/neuron_soma.h index 1c665c489..066fc8f81 100644 --- a/src/neuroscience/neuron_soma.h +++ b/src/neuroscience/neuron_soma.h @@ -28,7 +28,7 @@ namespace neuroscience { class NeuriteElement; class NeuronSoma : public Cell, public NeuronOrNeurite { - BDM_AGENT_HEADER(NeuronSoma, Cell, 1); + BDM_AGENT_HEADER(NeuronSoma, Cell); public: NeuronSoma(); diff --git a/src/neuroscience/param.h b/src/neuroscience/param.h index 29517da96..72966c24e 100644 --- a/src/neuroscience/param.h +++ b/src/neuroscience/param.h @@ -23,7 +23,6 @@ #include #include "core/param/param_group.h" #include "core/real_t.h" -#include "core/util/root.h" #include "cpptoml/cpptoml.h" namespace bdm { diff --git a/test/unit/core/agent/agent_pointer_test.cc b/test/unit/core/agent/agent_pointer_test.cc index 4bb772d5b..68e09be17 100644 --- a/test/unit/core/agent/agent_pointer_test.cc +++ b/test/unit/core/agent/agent_pointer_test.cc @@ -14,7 +14,7 @@ #include "unit/core/agent/agent_pointer_test.h" #include "core/randomized_rm.h" // for bdm::Ubrng -#include "unit/test_util/io_test.h" +#include "unit/test_util/test_util.h" namespace bdm { namespace agent_pointer_test_internal { @@ -123,7 +123,7 @@ void RunSortTest(Simulation* sim, AgentPointerMode mode) { } auto* random = Simulation::GetActive()->GetRandom(); - std::shuffle(ap_vector.begin(), ap_vector.end(), Ubrng(random)); + std::shuffle(ap_vector.begin(), ap_vector.end(), *random); std::sort(ap_vector.begin(), ap_vector.end()); for (uint64_t i = 0u; i < ap_vector.size(); ++i) { @@ -223,29 +223,5 @@ TEST(IsAgentPtrTest, All) { "AgentPointer is an AgentPointer"); } -#ifdef USE_DICT - -TEST_F(IOTest, AgentPointerIndirect) { - Simulation simulation(TEST_NAME); - RunIOTest(&simulation, AgentPointerMode::kIndirect); -} - -TEST_F(IOTest, AgentPointerDirect) { - Simulation simulation(TEST_NAME); - RunIOTest(&simulation, AgentPointerMode::kDirect); -} - -TEST_F(IOTest, AgentPointerNullptrIndirect) { - Simulation simulation(TEST_NAME); - IOTestAgentPointerNullptr(AgentPointerMode::kIndirect); -} - -TEST_F(IOTest, AgentPointerNullptrDirect) { - Simulation simulation(TEST_NAME); - IOTestAgentPointerNullptr(AgentPointerMode::kDirect); -} - -#endif // USE_DICT - } // namespace agent_pointer_test_internal } // namespace bdm diff --git a/test/unit/core/agent/agent_pointer_test.h b/test/unit/core/agent/agent_pointer_test.h index bd6c279b8..d082f72f4 100644 --- a/test/unit/core/agent/agent_pointer_test.h +++ b/test/unit/core/agent/agent_pointer_test.h @@ -20,48 +20,11 @@ #include "core/agent/agent_pointer.h" #include "core/resource_manager.h" #include "core/simulation.h" -#include "unit/test_util/io_test.h" #include "unit/test_util/test_agent.h" namespace bdm { namespace agent_pointer_test_internal { -inline void RunIOTest(Simulation* sim, AgentPointerMode mode) { - auto prev_mode = gAgentPointerMode; - gAgentPointerMode = mode; - - auto* rm = sim->GetResourceManager(); - rm->AddAgent(new TestAgent(123)); - TestAgent* so2 = new TestAgent(456); - rm->AddAgent(so2); - - AgentPointer agent_ptr(so2->GetUid()); - AgentPointer* restored; - - BackupAndRestore(agent_ptr, &restored); - - EXPECT_TRUE(*restored != nullptr); - EXPECT_EQ(456, (*restored)->GetData()); - - // restore gAgentPointerMode - gAgentPointerMode = prev_mode; -} - -inline void IOTestAgentPointerNullptr(AgentPointerMode mode) { - auto prev_mode = gAgentPointerMode; - gAgentPointerMode = mode; - - AgentPointer null_agent_pointer; - AgentPointer* restored = nullptr; - - BackupAndRestore(null_agent_pointer, &restored); - - EXPECT_TRUE(*restored == nullptr); - - // restore gAgentPointerMode - gAgentPointerMode = prev_mode; -} - } // namespace agent_pointer_test_internal } // namespace bdm diff --git a/test/unit/core/agent/agent_test.h b/test/unit/core/agent/agent_test.h index 3d3986236..5a246d252 100644 --- a/test/unit/core/agent/agent_test.h +++ b/test/unit/core/agent/agent_test.h @@ -25,7 +25,7 @@ namespace bdm { namespace agent_test_internal { struct Growth : public Behavior { - BDM_BEHAVIOR_HEADER(Growth, Behavior, 1); + BDM_BEHAVIOR_HEADER(Growth, Behavior); real_t growth_rate_ = 0.5; @@ -50,7 +50,7 @@ struct Growth : public Behavior { }; struct Movement : public Behavior { - BDM_BEHAVIOR_HEADER(Movement, Behavior, 1); + BDM_BEHAVIOR_HEADER(Movement, Behavior); Real3 velocity_; Movement() : velocity_({{0, 0, 0}}) { @@ -80,7 +80,7 @@ struct Movement : public Behavior { /// This behavior removes itself the first time it is executed struct Removal : public Behavior { - BDM_BEHAVIOR_HEADER(Removal, Movement, 1); + BDM_BEHAVIOR_HEADER(Removal, Movement); Removal() = default; virtual ~Removal() = default; @@ -92,7 +92,7 @@ struct Removal : public Behavior { // ----------------------------------------------------------------------------- struct CaptureStaticness : public Behavior { - BDM_BEHAVIOR_HEADER(CaptureStaticness, Behavior, 1); + BDM_BEHAVIOR_HEADER(CaptureStaticness, Behavior); CaptureStaticness() = default; CaptureStaticness(std::unordered_map* static_agents_map) @@ -108,10 +108,6 @@ struct CaptureStaticness : public Behavior { std::unordered_map* static_agents_map_; }; -#ifdef __ROOTCLING__ -static AgentPointer dummy_ptr; -#endif - } // namespace bdm #endif // UNIT_CORE_AGENT_AGENT_TEST_H_ diff --git a/test/unit/core/agent/agent_uid_generator_test.cc b/test/unit/core/agent/agent_uid_generator_test.cc index b044cbd2c..7476f53b8 100644 --- a/test/unit/core/agent/agent_uid_generator_test.cc +++ b/test/unit/core/agent/agent_uid_generator_test.cc @@ -16,7 +16,6 @@ #include #include "core/resource_manager.h" #include "core/simulation.h" -#include "unit/test_util/io_test.h" #include "unit/test_util/test_agent.h" namespace bdm { @@ -43,53 +42,4 @@ TEST(AgentUidGeneratorTest, NormalAndDefragmentationMode) { } } -#ifdef USE_DICT -TEST_F(IOTest, AgentUidGenerator) { - AgentUidGenerator test; - test.GenerateUid(); - test.GenerateUid(); - test.GenerateUid(); - - AgentUidGenerator* restored = nullptr; - - BackupAndRestore(test, &restored); - - EXPECT_EQ(restored->GetHighestIndex(), 3u); - EXPECT_EQ(restored->GenerateUid(), AgentUid(3u)); - - delete restored; -} - -TEST_F(IOTest, AgentUidGeneratorWithReuse) { - AgentUidGenerator generator; - - // Create num threads agent uids - auto* tinfo = ThreadInfo::GetInstance(); - for (int i = 0; i < 2 * tinfo->GetMaxThreads(); ++i) { - EXPECT_EQ(AgentUid(i), generator.GenerateUid()); - } - - // Mark for reuse half of the uids -#pragma omp parallel for schedule(static, 1) - for (int i = 0; i < tinfo->GetMaxThreads(); ++i) { - generator.ReuseAgentUid(AgentUid(i)); - } - - AgentUidGenerator* restored = nullptr; - - BackupAndRestore(generator, &restored); - - EXPECT_EQ(restored->GetHighestIndex(), 2u * tinfo->GetMaxThreads()); - - // Generate uids using the indices marked for reuse -#pragma omp parallel for schedule(static, 1) - for (int i = 0; i < tinfo->GetMaxThreads(); ++i) { - EXPECT_EQ(AgentUid(i, 1), restored->GenerateUid()); - } - - delete restored; -} - -#endif // USE_DICT - } // namespace bdm diff --git a/test/unit/core/agent/agent_uid_test.cc b/test/unit/core/agent/agent_uid_test.cc index 11c791738..428a2b3c5 100644 --- a/test/unit/core/agent/agent_uid_test.cc +++ b/test/unit/core/agent/agent_uid_test.cc @@ -14,7 +14,6 @@ #include "core/agent/agent_uid.h" #include -#include "unit/test_util/io_test.h" namespace bdm { @@ -98,21 +97,8 @@ TEST(AgentUidTest, uint64_tOperator) { TEST(AgentUidTest, uint64_tOperator2) { AgentUid uid(123, 2); uint64_t idx = uid; - EXPECT_EQ(idx, 8589934715u); // (2 << 32) | 123u); + constexpr uint64_t expected_uid = (uint64_t{2} << 32) | 123u; + EXPECT_EQ(idx, expected_uid); } -#ifdef USE_DICT -TEST_F(IOTest, AgentUid) { - AgentUid test{123u, 456u}; - AgentUid* restored = nullptr; - - BackupAndRestore(test, &restored); - - EXPECT_EQ(restored->GetIndex(), 123u); - EXPECT_EQ(restored->GetReused(), 456u); - - delete restored; -} -#endif // USE_DICT - } // namespace bdm diff --git a/test/unit/core/agent/cell_test.cc b/test/unit/core/agent/cell_test.cc index e11b8134e..7f185d47a 100644 --- a/test/unit/core/agent/cell_test.cc +++ b/test/unit/core/agent/cell_test.cc @@ -165,9 +165,5 @@ TEST(CellTest, DivideVolumeRatioAxis) { EXPECT_NEAR(cell.captured_theta_, 0.72664234068172562, kEpsilon); } -#ifdef USE_DICT -TEST(CellTest, IO) { RunIOTest(); } -#endif // USE_DICT - } // namespace cell_test_internal } // namespace bdm diff --git a/test/unit/core/agent/cell_test.h b/test/unit/core/agent/cell_test.h index 272abfcb4..de2d197d8 100644 --- a/test/unit/core/agent/cell_test.h +++ b/test/unit/core/agent/cell_test.h @@ -19,19 +19,16 @@ #include "core/agent/cell.h" #include "core/agent/cell_division_event.h" -#include "core/util/io.h" #include "gtest/gtest.h" #include "unit/core/agent/agent_test.h" #include "unit/test_util/test_util.h" -#define ROOTFILE "bdmFile.root" - namespace bdm { namespace cell_test_internal { /// Class used to get access to protected members class TestCell : public Cell { - BDM_AGENT_HEADER(TestCell, Cell, 1); + BDM_AGENT_HEADER(TestCell, Cell); public: TestCell() = default; @@ -70,68 +67,6 @@ class TestCell : public Cell { } // namespace cell_test_internal -namespace cell_test_internal { - -inline void RunIOTest() { - Simulation simulation("CellTest-RunIOTest"); - - using Growth = agent_test_internal::Growth; - using Movement = agent_test_internal::Movement; - remove(ROOTFILE); - - TestCell cell; - cell.SetPosition({5, 6, 7}); - cell.SetTractorForce({7, 4, 1}); - cell.SetDiameter(12); - cell.UpdateVolume(); - cell.SetAdherence(1.1); - cell.SetMass(5); - cell.AddBehavior(new Growth()); - cell.AddBehavior(new Movement({1, 2, 3})); - cell.SetBoxIdx(123); - - // write to root file - WritePersistentObject(ROOTFILE, "cell", cell, "new"); - - // read back - TestCell* restored_cell = nullptr; - GetPersistentObject(ROOTFILE, "cell", restored_cell); - - // validate - const real_t kEpsilon = abs_error::value; - EXPECT_NEAR(5, restored_cell->GetPosition()[0], kEpsilon); - EXPECT_NEAR(6, restored_cell->GetPosition()[1], kEpsilon); - EXPECT_NEAR(7, restored_cell->GetPosition()[2], kEpsilon); - - EXPECT_NEAR(7, restored_cell->GetTractorForce()[0], kEpsilon); - EXPECT_NEAR(4, restored_cell->GetTractorForce()[1], kEpsilon); - EXPECT_NEAR(1, restored_cell->GetTractorForce()[2], kEpsilon); - - EXPECT_NEAR(12, restored_cell->GetDiameter(), kEpsilon); - // differs slightly from the value in branch validation due to more precise - // value of PI - EXPECT_NEAR(cell.GetVolume(), restored_cell->GetVolume(), kEpsilon); - EXPECT_NEAR(1.1, restored_cell->GetAdherence(), kEpsilon); - EXPECT_NEAR(5, restored_cell->GetMass(), kEpsilon); - - EXPECT_EQ(2u, restored_cell->GetAllBehaviors().size()); - EXPECT_TRUE(dynamic_cast(restored_cell->GetAllBehaviors()[0]) != - nullptr); - EXPECT_NEAR( - 0.5, - dynamic_cast(restored_cell->GetAllBehaviors()[0])->growth_rate_, - kEpsilon); - EXPECT_TRUE(dynamic_cast(restored_cell->GetAllBehaviors()[1]) != - nullptr); - - EXPECT_EQ(123u, restored_cell->GetBoxIdx()); - - delete restored_cell; - // delete root file - remove(ROOTFILE); -} - -} // namespace cell_test_internal } // namespace bdm #endif // UNIT_CORE_AGENT_CELL_TEST_H_ diff --git a/test/unit/core/analysis/reduce_test.cc b/test/unit/core/analysis/reduce_test.cc index b06020827..6e0636751 100644 --- a/test/unit/core/analysis/reduce_test.cc +++ b/test/unit/core/analysis/reduce_test.cc @@ -17,19 +17,21 @@ #include "core/resource_manager.h" #include "core/scheduler.h" #include "core/simulation.h" -#include "unit/test_util/io_test.h" #include "unit/test_util/test_agent.h" #include "unit/test_util/test_util.h" namespace bdm { namespace experimental { +constexpr uint64_t kAgentCount = 2000; +constexpr uint64_t kAgentDataSum = kAgentCount * (kAgentCount - 1) / 2; + // ----------------------------------------------------------------------------- TEST(Reduce, Reduce) { Simulation sim(TEST_NAME); auto* rm = sim.GetResourceManager(); - for (uint64_t i = 0; i < 2000; ++i) { + for (uint64_t i = 0; i < kAgentCount; ++i) { auto* a = new TestAgent(); a->SetData(i); rm->AddAgent(a); @@ -40,7 +42,7 @@ TEST(Reduce, Reduce) { }); SumReduction combine_tl_results; auto result = Reduce(&sim, sum_data, combine_tl_results); - EXPECT_EQ(1999000u, result); + EXPECT_EQ(kAgentDataSum, result); } // ----------------------------------------------------------------------------- @@ -48,7 +50,7 @@ TEST(Reduce, GenericReducer) { Simulation sim(TEST_NAME); auto* rm = sim.GetResourceManager(); - for (uint64_t i = 0; i < 2000; ++i) { + for (uint64_t i = 0; i < kAgentCount; ++i) { auto* a = new TestAgent(); a->SetData(i); rm->AddAgent(a); @@ -70,7 +72,7 @@ TEST(Reduce, GenericReducer) { GenericReducer reducer(sum_data, combine_tl_results); rm->ForEachAgentParallel(reducer); auto result = reducer.GetResult(); - EXPECT_EQ(1999000u, result); + EXPECT_EQ(kAgentDataSum, result); } // with filter, without post processing @@ -181,69 +183,5 @@ TEST(Reduce, Counter) { } } -#ifdef USE_DICT -// ----------------------------------------------------------------------------- -TEST_F(IOTest, GenericReducer) { - Simulation sim(TEST_NAME); - auto* rm = sim.GetResourceManager(); - - for (uint64_t i = 0; i < 2000; ++i) { - auto* a = new TestAgent(); - a->SetData(i); - rm->AddAgent(a); - } - - auto sum_data = [](Agent* agent, uint64_t* tl_result) { - *tl_result += bdm_static_cast(agent)->GetData(); - }; - auto combine_tl_results = [](const SharedData& tl_results) { - uint64_t result = 0; - for (auto& el : tl_results) { - result += el; - } - return result; - }; - auto post_process = [](uint64_t result) { return result / 2; }; - auto filter = [](Agent* a) { - return bdm_static_cast(a)->GetData() <= 1000; - }; - GenericReducer reducer(sum_data, combine_tl_results, filter, - post_process); - - GenericReducer* restored; - BackupAndRestore(reducer, &restored); - - rm->ForEachAgentParallel(*restored); - auto result = restored->GetResult(); - EXPECT_EQ(250250u, result); -} - -// ----------------------------------------------------------------------------- -TEST_F(IOTest, Counter) { - Simulation sim(TEST_NAME); - auto* rm = sim.GetResourceManager(); - - for (uint64_t i = 0; i < 2000; ++i) { - auto* a = new TestAgent(); - a->SetData(i); - rm->AddAgent(a); - } - - auto data_lt_1000 = [](Agent* agent) { - return bdm_static_cast(agent)->GetData() < 1000; - }; - auto post_process = [](uint64_t result) { return result / 2; }; - Counter<> counter(data_lt_1000, post_process); - - Counter<>* restored; - BackupAndRestore(counter, &restored); - - rm->ForEachAgentParallel(*restored); - auto result = restored->GetResult(); - EXPECT_EQ(500u, result); -} - -#endif // USE_DICT - } // namespace experimental } // namespace bdm diff --git a/test/unit/core/behavior/stateless_behavior_test.cc b/test/unit/core/behavior/stateless_behavior_test.cc index 9cd4df4a3..b96fe535b 100644 --- a/test/unit/core/behavior/stateless_behavior_test.cc +++ b/test/unit/core/behavior/stateless_behavior_test.cc @@ -15,7 +15,6 @@ #include "core/behavior/stateless_behavior.h" #include #include "core/agent/cell_division_event.h" -#include "unit/test_util/io_test.h" #include "unit/test_util/test_agent.h" #include "unit/test_util/test_util.h" @@ -109,22 +108,4 @@ TEST(StatelessBehavior, Event) { } } -#ifdef USE_DICT -// ----------------------------------------------------------------------------- -TEST_F(IOTest, StatelessBehavior) { - Simulation sim(TEST_NAME); - - StatelessBehavior b( - [](Agent* a) { bdm_static_cast(a)->SetData(123); }); - - StatelessBehavior* restored; - BackupAndRestore(b, &restored); - - TestAgent a; - restored->Run(&a); - EXPECT_EQ(123, a.GetData()); -} - -#endif // USE_DICT - } // namespace bdm diff --git a/test/unit/test_util/test_agent.h b/test/unit/test_util/test_agent.h index 10e854e2c..6ddaee32d 100644 --- a/test/unit/test_util/test_agent.h +++ b/test/unit/test_util/test_agent.h @@ -22,7 +22,7 @@ namespace bdm { class TestAgent : public Agent { - BDM_AGENT_HEADER(TestAgent, Agent, 1); + BDM_AGENT_HEADER(TestAgent, Agent); public: TestAgent() = default;