Reuse trajectories for multiple experiments

src/dynamics/CMakeLists.txt

@@ -1,4 +1,4 @@
-add_library(dynamics STATIC
-  base.h
+add_library(dynamics OBJECT
+  base.cpp base.h
   decoupled.cpp decoupled.h
 )

src/dynamics/base.cpp

@@ -0,0 +1,6 @@
+#include "base.h"
+#include "../utils/registry.h"
+
+std::unique_ptr<Dynamics> Dynamics::createFromInput(const std::string &input) {
+  return Registry<Dynamics>::create(input);
+}

src/dynamics/base.h

@@ -22,6 +22,8 @@ public:
   [[nodiscard]] virtual double getInitOmega() const = 0;
   [[nodiscard]] virtual std::unique_ptr<Dynamics> clone() const = 0;
   [[nodiscard]] virtual std::string toString() const = 0;
+
+  static std::unique_ptr<Dynamics> createFromInput(const std::string &input);
 };
 
 #endif // RWSIM_DYNAMICSBASE_H

src/main.cpp

@@ -1,4 +1,5 @@
 
+#include "dynamics/base.h"
 #include "dynamics/decoupled.h"
 #include "experiments/spectrum.h"
 #include "experiments/ste.h"
@@ -15,7 +16,7 @@ int main(const int argc, char *argv[]) {
   Arguments args;
   try {
     args = parse_args(argc, argv);
-  } catch (std::invalid_argument &error) {
+  } catch (std::exception &error) {
     std::cerr << error.what() << std::endl;
     return 1;
   }
@@ -41,17 +42,28 @@ int main(const int argc, char *argv[]) {
     rng.seed(rd());
   }
 
-  DecoupledDynamics dynamics(
-      motions::BaseMotion::createFromInput(args.motion_type),
-      times::BaseDistribution::createFromInput(args.distribution_type));
+  std::unique_ptr<Dynamics> dynamics;
+  if (!args.dynamics_type.empty()) {
+    dynamics = Dynamics::createFromInput(args.dynamics_type);
+  } else {
+    dynamics = std::make_unique<DecoupledDynamics>(
+        motions::BaseMotion::createFromInput(args.motion_type),
+        times::BaseDistribution::createFromInput(args.distribution_type));
+  }
 
+  SpectrumExperiment spectrum_exp;
+  STEExperiment ste_exp;
+
+  std::vector<Experiment *> experiments;
   if (args.spectrum) {
-    SpectrumExperiment experiment;
-    run_simulation(experiment, parameter, args.optional, dynamics, rng);
+    experiments.push_back(&spectrum_exp);
   }
   if (args.ste) {
-    STEExperiment experiment;
-    run_simulation(experiment, parameter, args.optional, dynamics, rng);
+    experiments.push_back(&ste_exp);
+  }
+
+  if (!experiments.empty()) {
+    run_simulation(experiments, parameter, args.optional, *dynamics, rng);
   }
 
   return 0;

src/sims.cpp

@@ -2,18 +2,24 @@
 #include "utils/functions.h"
 #include "utils/io.h"
 
+#include <algorithm>
 #include <chrono>
 #include <iostream>
 #include <omp.h>
 
-void run_simulation(Experiment &experiment,
+void run_simulation(std::vector<Experiment *> experiments,
                     const std::unordered_map<std::string, double> &parameter,
                     const std::unordered_map<std::string, double> &optional,
                     Dynamics &dynamics, std::mt19937_64 &rng) {
   const int num_walker = static_cast<int>(parameter.at("num_walker"));
 
   dynamics.setParameters(parameter);
-  experiment.setup(parameter, optional);
+
+  double tmax = 0;
+  for (auto *exp : experiments) {
+    exp->setup(parameter, optional);
+    tmax = std::max(tmax, exp->tmax());
+  }
 
   const auto start = printStart(optional);
 
@@ -26,12 +32,15 @@ void run_simulation(Experiment &experiment,
     thread_rngs.emplace_back(rng());
   }
 
-  // Create per-thread clones of dynamics and experiment
+  // Create per-thread clones of dynamics and experiments
   std::vector<std::unique_ptr<Dynamics>> thread_dynamics;
-  std::vector<std::unique_ptr<Experiment>> thread_experiments;
+  std::vector<std::vector<std::unique_ptr<Experiment>>> thread_experiments(
+      num_threads);
   for (int i = 0; i < num_threads; i++) {
     thread_dynamics.push_back(dynamics.clone());
-    thread_experiments.push_back(experiment.clone());
+    for (auto *exp : experiments) {
+      thread_experiments[i].push_back(exp->clone());
+    }
   }
 
   int steps_done = 0;
@@ -42,13 +51,16 @@ void run_simulation(Experiment &experiment,
     const int tid = omp_get_thread_num();
     auto &local_rng = thread_rngs[tid];
     auto &local_dynamics = *thread_dynamics[tid];
-    auto &local_experiment = *thread_experiments[tid];
+    auto &local_experiments = thread_experiments[tid];
 
 #pragma omp for schedule(static)
     for (int mol_i = 0; mol_i < num_walker; mol_i++) {
-      auto traj = make_trajectory(local_dynamics, experiment.tmax(), local_rng);
-      local_experiment.accumulate(traj, local_dynamics.getInitOmega(),
-                                  num_walker);
+      auto traj = make_trajectory(local_dynamics, tmax, local_rng);
+      const double init_omega = local_dynamics.getInitOmega();
+
+      for (auto &exp : local_experiments) {
+        exp->accumulate(traj, init_omega, num_walker);
+      }
 
       if (tid == 0) {
 #pragma omp atomic
@@ -62,13 +74,15 @@ void run_simulation(Experiment &experiment,
     }
   }
 
-  // Merge per-thread results
-  for (int i = 0; i < num_threads; i++) {
-    experiment.merge(*thread_experiments[i]);
+  // Merge per-thread results and save
+  const auto directory = make_directory(dynamics.toString());
+  for (size_t e = 0; e < experiments.size(); e++) {
+    for (int i = 0; i < num_threads; i++) {
+      experiments[e]->merge(*thread_experiments[i][e]);
+    }
+    experiments[e]->save(directory);
   }
 
-  const auto directory = make_directory(dynamics.toString());
-  experiment.save(directory);
   printEnd(start);
 }
 

src/sims.h

@@ -8,8 +8,9 @@
 #include <random>
 #include <string>
 #include <unordered_map>
+#include <vector>
 
-void run_simulation(Experiment &experiment,
+void run_simulation(std::vector<Experiment *> experiments,
                     const std::unordered_map<std::string, double> &parameter,
                     const std::unordered_map<std::string, double> &optional,
                     Dynamics &dynamics, std::mt19937_64 &rng);

src/utils/io.cpp

@@ -41,6 +41,9 @@ Arguments parse_args(const int argc, char *argv[]) {
   // convert to vector to use iterator for loop
   const std::vector<std::string> input_args(argv + 1, argv + argc);
 
+  // collect positional arguments separately
+  std::vector<std::string> positional;
+
   for (auto it = input_args.begin(); it != input_args.end(); ++it) {
     // check for optional parameter
     if (it->at(0) == '-') {
@@ -62,29 +65,23 @@ Arguments parse_args(const int argc, char *argv[]) {
       continue;
     }
 
-    // Two positional parameters are defined: 1. Location of config file; 2.
-    // Name of motion model
-    if (args.parameter_file.empty()) {
-      args.parameter_file = *it;
-      continue;
-    }
-
-    if (args.motion_type.empty()) {
-      args.motion_type = *it;
-      continue;
-    }
-
-    if (args.distribution_type.empty()) {
-      args.distribution_type = *it;
-      continue;
-    }
-
-    throw std::invalid_argument("too many positional arguments");
+    positional.push_back(*it);
   }
 
-  if (args.parameter_file.empty() || args.motion_type.empty() ||
-      args.distribution_type.empty()) {
-    throw std::invalid_argument("Missing argument");
+  if (positional.size() == 2) {
+    // config_file DynamicsModel
+    args.parameter_file = positional[0];
+    args.dynamics_type = positional[1];
+  } else if (positional.size() == 3) {
+    // config_file MotionModel DistributionModel
+    args.parameter_file = positional[0];
+    args.motion_type = positional[1];
+    args.distribution_type = positional[2];
+  } else {
+    throw std::invalid_argument(
+        "Expected 2 positional args (config dynamics) or "
+        "3 positional args (config motion distribution), got " +
+        std::to_string(positional.size()));
   }
 
   return args;

src/utils/io.h

@@ -13,6 +13,7 @@ struct Arguments {
   bool spectrum = false;
   std::string motion_type{};
   std::string distribution_type{};
+  std::string dynamics_type{};
   std::unordered_map<std::string, double> optional;
 };