Create Experiment class

src/CMakeLists.txt

@@ -2,15 +2,16 @@
 add_subdirectory(times)
 add_subdirectory(motions)
 add_subdirectory(utils)
+add_subdirectory(experiments)
 
 add_library(simulation STATIC sims.cpp sims.h)
-target_link_libraries(simulation PRIVATE utils)
+target_link_libraries(simulation PRIVATE utils experiments)
 
 add_executable(
         rwsim
         main.cpp
 )
 
-target_link_libraries(rwsim PUBLIC RWMotion RWTime utils simulation)
+target_link_libraries(rwsim PUBLIC RWMotion RWTime utils experiments simulation)
 
 target_compile_options(rwsim PUBLIC -Werror -Wall -Wextra -Wconversion -O2)

src/experiments/CMakeLists.txt

@@ -0,0 +1,10 @@
+add_library(
+        experiments STATIC
+        base.h
+        spectrum.cpp
+        spectrum.h
+        ste.cpp
+        ste.h
+)
+
+target_link_libraries(experiments PRIVATE utils)

src/experiments/base.h

@@ -0,0 +1,31 @@
+#ifndef RWSIM_EXPERIMENTBASE_H
+#define RWSIM_EXPERIMENTBASE_H
+
+#include "../motions/base.h"
+#include "../times/base.h"
+
+#include <memory>
+#include <unordered_map>
+#include <string>
+#include <vector>
+
+struct Trajectory {
+    std::vector<double> time;
+    std::vector<double> phase;
+    std::vector<double> omega;
+};
+
+class Experiment {
+public:
+    virtual ~Experiment() = default;
+
+    virtual void setup(const std::unordered_map<std::string, double>& parameter,
+                       const std::unordered_map<std::string, double>& optional) = 0;
+    [[nodiscard]] virtual double tmax() const = 0;
+    virtual void accumulate(const Trajectory& traj, double init_omega, int num_walkers) = 0;
+    virtual void save(const motions::BaseMotion& motion, const times::BaseDistribution& dist) = 0;
+    [[nodiscard]] virtual std::unique_ptr<Experiment> clone() const = 0;
+    virtual void merge(const Experiment& other) = 0;
+};
+
+#endif //RWSIM_EXPERIMENTBASE_H

src/experiments/spectrum.cpp

@@ -0,0 +1,65 @@
+#include "spectrum.h"
+#include "../utils/functions.h"
+#include "../utils/ranges.h"
+#include "../utils/io.h"
+
+#include <algorithm>
+#include <cmath>
+
+void SpectrumExperiment::setup(
+    const std::unordered_map<std::string, double>& parameter,
+    const std::unordered_map<std::string, double>& optional
+    ) {
+    m_parameter = parameter;
+    m_optional = optional;
+    m_num_acq = static_cast<int>(parameter.at("num_acq"));
+    m_t_fid = arange(m_num_acq, parameter.at("dwell_time"));
+    m_echo_times = linspace(parameter.at("techo_start"), parameter.at("techo_stop"), static_cast<int>(parameter.at("techo_steps")));
+
+    m_fid_dict.clear();
+    for (auto t_echo_i: m_echo_times) {
+        m_fid_dict[t_echo_i] = std::vector<double>(m_num_acq, 0.);
+    }
+
+    m_tmax = *std::max_element(m_echo_times.begin(), m_echo_times.end()) * 2 + m_t_fid.back();
+}
+
+double SpectrumExperiment::tmax() const {
+    return m_tmax;
+}
+
+void SpectrumExperiment::accumulate(const Trajectory& traj, double, int num_walkers) {
+    for (auto& [t_echo_j, fid_j] : m_fid_dict) {
+        int current_pos = nearest_index(traj.time, t_echo_j, 0);
+        const double phase_techo = lerp(traj.time, traj.phase, t_echo_j, current_pos);
+
+        for (int acq_idx = 0; acq_idx < m_num_acq; acq_idx++) {
+            const double real_time = m_t_fid[acq_idx] + 2 * t_echo_j;
+
+            current_pos = nearest_index(traj.time, real_time, current_pos);
+            const double phase_acq = lerp(traj.time, traj.phase, real_time, current_pos);
+
+            fid_j[acq_idx] += std::cos(phase_acq - 2 * phase_techo) / num_walkers;
+        }
+    }
+}
+
+void SpectrumExperiment::save(const motions::BaseMotion& motion, const times::BaseDistribution& dist) {
+    const auto path = make_directory(motion, dist);
+    save_parameter_to_file(std::string("timesignal"), path, m_parameter, m_optional);
+    save_data_to_file(std::string("timesignal"), path, m_t_fid, m_fid_dict, m_optional);
+}
+
+std::unique_ptr<Experiment> SpectrumExperiment::clone() const {
+    return std::make_unique<SpectrumExperiment>(*this);
+}
+
+void SpectrumExperiment::merge(const Experiment& other) {
+    const auto& o = dynamic_cast<const SpectrumExperiment&>(other);
+    for (auto& [t_echo, fid] : m_fid_dict) {
+        const auto& other_fid = o.m_fid_dict.at(t_echo);
+        for (size_t i = 0; i < fid.size(); i++) {
+            fid[i] += other_fid[i];
+        }
+    }
+}

src/experiments/spectrum.h

@@ -0,0 +1,30 @@
+#ifndef RWSIM_SPECTRUM_H
+#define RWSIM_SPECTRUM_H
+
+#include "base.h"
+
+#include <map>
+#include <vector>
+#include <unordered_map>
+
+class SpectrumExperiment final : public Experiment {
+public:
+    void setup(const std::unordered_map<std::string, double>& parameter,
+               const std::unordered_map<std::string, double>& optional) override;
+    [[nodiscard]] double tmax() const override;
+    void accumulate(const Trajectory& traj, double init_omega, int num_walkers) override;
+    void save(const motions::BaseMotion& motion, const times::BaseDistribution& dist) override;
+    [[nodiscard]] std::unique_ptr<Experiment> clone() const override;
+    void merge(const Experiment& other) override;
+
+private:
+    int m_num_acq{0};
+    std::vector<double> m_t_fid;
+    std::vector<double> m_echo_times;
+    std::map<double, std::vector<double>> m_fid_dict;
+    std::unordered_map<std::string, double> m_parameter;
+    std::unordered_map<std::string, double> m_optional;
+    double m_tmax{0};
+};
+
+#endif //RWSIM_SPECTRUM_H

src/experiments/ste.cpp

@@ -0,0 +1,99 @@
+#include "ste.h"
+#include "../utils/functions.h"
+#include "../utils/ranges.h"
+#include "../utils/io.h"
+
+#include <algorithm>
+#include <cmath>
+
+void STEExperiment::setup(
+    const std::unordered_map<std::string, double>& parameter,
+    const std::unordered_map<std::string, double>& optional
+    ) {
+    m_parameter = parameter;
+    m_optional = optional;
+    m_num_mix_times = static_cast<int>(parameter.at("tmix_steps"));
+    m_evolution_times = linspace(parameter.at("tevo_start"), parameter.at("tevo_stop"), static_cast<int>(parameter.at("tevo_steps")));
+    m_mixing_times = logspace(parameter.at("tmix_start"), parameter.at("tmix_stop"), m_num_mix_times);
+    m_tpulse4 = parameter.at("tpulse4");
+
+    m_cc_dict.clear();
+    m_ss_dict.clear();
+    for (auto t_evo_i: m_evolution_times) {
+        m_cc_dict[t_evo_i] = std::vector<double>(m_num_mix_times, 0.);
+        m_ss_dict[t_evo_i] = std::vector<double>(m_num_mix_times, 0.);
+    }
+    m_f2.assign(m_num_mix_times, 0.);
+
+    m_tmax = *std::max_element(m_evolution_times.begin(), m_evolution_times.end()) * 2
+             + *std::max_element(m_mixing_times.begin(), m_mixing_times.end())
+             + 2 * m_tpulse4;
+}
+
+double STEExperiment::tmax() const {
+    return m_tmax;
+}
+
+void STEExperiment::accumulate(const Trajectory& traj, double init_omega, int num_walkers) {
+    int f2_pos = 0;
+    for (int f2_idx = 0; f2_idx < m_num_mix_times; f2_idx++) {
+        const double t_mix_f2 = m_mixing_times[f2_idx];
+        f2_pos = nearest_index(traj.time, t_mix_f2, f2_pos);
+        m_f2[f2_idx] += traj.omega[f2_pos] * init_omega / num_walkers;
+    }
+
+    for (auto& [t_evo_j, _] : m_cc_dict) {
+        auto& cc_j = m_cc_dict[t_evo_j];
+        auto& ss_j = m_ss_dict[t_evo_j];
+
+        int current_pos = nearest_index(traj.time, t_evo_j, 0);
+        const double dephased = lerp(traj.time, traj.phase, t_evo_j, current_pos);
+        const double cc_tevo = std::cos(dephased);
+        const double ss_tevo = std::sin(dephased);
+
+        for (int mix_idx = 0; mix_idx < m_num_mix_times; mix_idx++) {
+            const double time_end_mix = m_mixing_times[mix_idx] + t_evo_j;
+            current_pos = nearest_index(traj.time, time_end_mix, current_pos);
+            const double phase_mix_end = lerp(traj.time, traj.phase, time_end_mix, current_pos);
+
+            const double time_pulse4 = time_end_mix + m_tpulse4;
+            current_pos = nearest_index(traj.time, time_pulse4, current_pos);
+            const double phase_4pulse = lerp(traj.time, traj.phase, time_pulse4, current_pos);
+
+            const double time_echo = time_pulse4 + m_tpulse4 + t_evo_j;
+            current_pos = nearest_index(traj.time, time_echo, current_pos);
+            double rephased = lerp(traj.time, traj.phase, time_echo, current_pos) + phase_mix_end - 2 * phase_4pulse;
+
+            cc_j[mix_idx] += cc_tevo * std::cos(rephased) / num_walkers;
+            ss_j[mix_idx] += ss_tevo * std::sin(rephased) / num_walkers;
+        }
+    }
+}
+
+void STEExperiment::save(const motions::BaseMotion& motion, const times::BaseDistribution& dist) {
+    const auto folders = make_directory(motion, dist);
+    save_parameter_to_file(std::string("ste"), folders, m_parameter, m_optional);
+    save_data_to_file(std::string("coscos"), folders, m_mixing_times, m_cc_dict, m_optional);
+    save_data_to_file(std::string("sinsin"), folders, m_mixing_times, m_ss_dict, m_optional);
+    save_data_to_file(std::string("f2"), folders, m_mixing_times, m_f2, m_optional);
+}
+
+std::unique_ptr<Experiment> STEExperiment::clone() const {
+    return std::make_unique<STEExperiment>(*this);
+}
+
+void STEExperiment::merge(const Experiment& other) {
+    const auto& o = dynamic_cast<const STEExperiment&>(other);
+    for (auto& [t_evo, cc] : m_cc_dict) {
+        const auto& other_cc = o.m_cc_dict.at(t_evo);
+        const auto& other_ss = o.m_ss_dict.at(t_evo);
+        auto& ss = m_ss_dict[t_evo];
+        for (size_t i = 0; i < cc.size(); i++) {
+            cc[i] += other_cc[i];
+            ss[i] += other_ss[i];
+        }
+    }
+    for (size_t i = 0; i < m_f2.size(); i++) {
+        m_f2[i] += o.m_f2[i];
+    }
+}

src/experiments/ste.h

@@ -0,0 +1,33 @@
+#ifndef RWSIM_STE_H
+#define RWSIM_STE_H
+
+#include "base.h"
+
+#include <map>
+#include <vector>
+#include <unordered_map>
+
+class STEExperiment final : public Experiment {
+public:
+    void setup(const std::unordered_map<std::string, double>& parameter,
+               const std::unordered_map<std::string, double>& optional) override;
+    [[nodiscard]] double tmax() const override;
+    void accumulate(const Trajectory& traj, double init_omega, int num_walkers) override;
+    void save(const motions::BaseMotion& motion, const times::BaseDistribution& dist) override;
+    [[nodiscard]] std::unique_ptr<Experiment> clone() const override;
+    void merge(const Experiment& other) override;
+
+private:
+    int m_num_mix_times{0};
+    double m_tpulse4{0};
+    std::vector<double> m_evolution_times;
+    std::vector<double> m_mixing_times;
+    std::map<double, std::vector<double>> m_cc_dict;
+    std::map<double, std::vector<double>> m_ss_dict;
+    std::vector<double> m_f2;
+    std::unordered_map<std::string, double> m_parameter;
+    std::unordered_map<std::string, double> m_optional;
+    double m_tmax{0};
+};
+
+#endif //RWSIM_STE_H

src/main.cpp

@@ -1,5 +1,7 @@
 
 #include "sims.h"
+#include "experiments/spectrum.h"
+#include "experiments/ste.h"
 #include "utils/io.h"
 #include "motions/base.h"
 #include "times/base.h"
@@ -43,11 +45,14 @@ int main (const int argc, char *argv[]) {
 
     auto motion = motions::BaseMotion::createFromInput(args.motion_type);
     auto dist = times::BaseDistribution::createFromInput(args.distribution_type);
+
     if (args.spectrum) {
-       run_spectrum(parameter, args.optional, *motion, *dist, rng);
+        SpectrumExperiment experiment;
+        run_simulation(experiment, parameter, args.optional, *motion, *dist, rng);
     }
     if (args.ste) {
-        run_ste(parameter, args.optional, *motion, *dist, rng);
+        STEExperiment experiment;
+        run_simulation(experiment, parameter, args.optional, *motion, *dist, rng);
     }
 
     return 0;

src/sims.cpp

@@ -1,169 +1,34 @@
 #include "sims.h"
-#include "times/base.h"
 #include "utils/functions.h"
-#include "utils/ranges.h"
-#include "utils/io.h"
 
 #include <iostream>
-#include <algorithm>
-#include <unordered_map>
-#include <map>
-#include <string>
-#include <vector>
-#include <cmath>
 #include <chrono>
 
 
-
+void run_simulation(
-void run_spectrum(
+    Experiment& experiment,
     std::unordered_map<std::string, double>& parameter,
-    std::unordered_map<std::string, double> optional,
+    std::unordered_map<std::string, double>& optional,
     motions::BaseMotion& motion,
     times::BaseDistribution& dist,
     std::mt19937_64& rng
     ) {
     const int num_walker = static_cast<int>(parameter["num_walker"]);
 
-    // time axis for all time signals
-    const int num_acq = static_cast<int>(parameter["num_acq"]);
-    const std::vector<double> t_fid = arange(num_acq, parameter["dwell_time"]);
-    const std::vector<double> echo_times = linspace(parameter["techo_start"], parameter["techo_stop"], static_cast<int>(parameter["techo_steps"]));
-
-    // make timesignal vectors and set them to zero
-    std::map<double, std::vector<double>> fid_dict;
-    for (auto t_echo_i: echo_times) {
-        fid_dict[t_echo_i] = std::vector<double>(num_acq);
-        std::fill(fid_dict[t_echo_i].begin(), fid_dict[t_echo_i].end(), 0.);
-    }
-
-    // calculate min length of a trajectory
-    const double tmax = *std::max_element(echo_times.begin(), echo_times.end()) * 2 + t_fid.back();
-
-    // set parameter in distribution and motion model
     dist.setParameters(parameter);
     motion.setParameters(parameter);
+    experiment.setup(parameter, optional);
 
     const auto start = printStart(optional);
     auto last_print_out = std::chrono::system_clock::now();
 
-    // let the walker walk
     for (int mol_i = 0; mol_i < num_walker; mol_i++) {
-        auto traj = make_trajectory(motion, dist, tmax, rng);
+        auto traj = make_trajectory(motion, dist, experiment.tmax(), rng);
-
+        experiment.accumulate(traj, motion.getInitOmega(), num_walker);
-        for (auto& [t_echo_j, fid_j] : fid_dict) {
-            // get phase at echo pulse
-            int current_pos = nearest_index(traj.time, t_echo_j, 0);
-            const double phase_techo = lerp(traj.time, traj.phase, t_echo_j, current_pos);
-
-
-            for (int acq_idx = 0; acq_idx < num_acq; acq_idx++) {
-                const double real_time = t_fid[acq_idx] + 2 * t_echo_j;
-
-                current_pos = nearest_index(traj.time, real_time, current_pos);
-                const double phase_acq = lerp(traj.time, traj.phase, real_time, current_pos);
-
-                fid_j[acq_idx] += std::cos(phase_acq - 2 * phase_techo) / num_walker;
-            }
-            last_print_out = printSteps(last_print_out, start, num_walker, mol_i);
-        }
-    }
-
-    // write fid to files
-    const auto path = make_directory(motion, dist);
-    save_parameter_to_file(std::string("timesignal"), path, parameter, optional);
-    save_data_to_file(std::string("timesignal"), path, t_fid, fid_dict, optional);
-
-    printEnd(start);
-}
-
-
-void run_ste(
-    std::unordered_map<std::string, double>& parameter,
-    std::unordered_map<std::string, double> optional,
-    motions::BaseMotion& motion,
-    times::BaseDistribution& dist,
-    std::mt19937_64& rng
-    ) {
-    const int num_walker = static_cast<int>(parameter[std::string("num_walker")]);
-
-    const int num_mix_times = static_cast<int>(parameter["tmix_steps"]);
-    const std::vector<double> evolution_times = linspace(parameter["tevo_start"], parameter["tevo_stop"], static_cast<int>(parameter["tevo_steps"]));
-    const std::vector<double> mixing_times = logspace(parameter["tmix_start"], parameter["tmix_stop"], num_mix_times);
-    const double tpulse4 = parameter["tpulse4"];
-
-    // make ste decay vectors and set them to zero
-    std::map<double, std::vector<double>> cc_dict;
-    std::map<double, std::vector<double>> ss_dict;
-    for (auto t_evo_i: evolution_times) {
-        cc_dict[t_evo_i] = std::vector<double>(num_mix_times);
-        ss_dict[t_evo_i] = std::vector<double>(num_mix_times);
-        std::fill(cc_dict[t_evo_i].begin(), cc_dict[t_evo_i].end(), 0.);
-        std::fill(ss_dict[t_evo_i].begin(), ss_dict[t_evo_i].end(), 0.);
-    }
-    std::vector<double> f2(num_mix_times);
-
-    // each trajectory must have a duration of at least tmax
-    const double tmax = *std::max_element(evolution_times.begin(), evolution_times.end()) * 2 + *std::max_element(mixing_times.begin(), mixing_times.end()) + 2*tpulse4;
-
-    // set parameter in distribution and motion model
-    dist.setParameters(parameter);
-    motion.setParameters(parameter);
-
-    const auto start = printStart(optional);
-    auto last_print_out = std::chrono::system_clock::now();
-
-    // let the walker walk
-    for (int mol_i = 0; mol_i < num_walker; mol_i++){
-        auto traj = make_trajectory(motion, dist, tmax, rng);
-
-        int f2_pos = 0;
-        for (int f2_idx=0; f2_idx < num_mix_times; f2_idx++) {
-            const double t_mix_f2 = mixing_times[f2_idx];
-            f2_pos = nearest_index(traj.time, t_mix_f2, f2_pos);
-            f2[f2_idx] += traj.omega[f2_pos] * motion.getInitOmega() / num_walker;
-        }
-
-
-        for (auto& [t_evo_j, _] : cc_dict) {
-            auto& cc_j = cc_dict[t_evo_j];
-            auto& ss_j = ss_dict[t_evo_j];
-
-            // get phase at beginning of mixing time
-            int current_pos = nearest_index(traj.time, t_evo_j, 0);
-            const double dephased = lerp(traj.time, traj.phase, t_evo_j, current_pos);
-            const double cc_tevo = std::cos(dephased);
-            const double ss_tevo = std::sin(dephased);
-
-            for (int mix_idx = 0; mix_idx < num_mix_times; mix_idx++) {
-                // get phase at end of mixing time
-                const double time_end_mix = mixing_times[mix_idx] + t_evo_j;
-                current_pos = nearest_index(traj.time, time_end_mix, current_pos);
-                const double phase_mix_end = lerp(traj.time, traj.phase, time_end_mix, current_pos);
-
-                // get phase at position of 4th pulse
-                const double time_pulse4 = time_end_mix + tpulse4;
-                current_pos = nearest_index(traj.time, time_pulse4, current_pos);
-                const double phase_4pulse = lerp(traj.time, traj.phase, time_pulse4, current_pos);
-
-                // get phase at echo position
-                const double time_echo = time_pulse4 + tpulse4 + t_evo_j;
-                current_pos = nearest_index(traj.time, time_echo, current_pos);
-                double rephased = lerp(traj.time, traj.phase, time_echo, current_pos) + phase_mix_end - 2*phase_4pulse;
-
-                cc_j[mix_idx] += cc_tevo * std::cos(rephased) / num_walker;
-                ss_j[mix_idx] += ss_tevo * std::sin(rephased) / num_walker;
-            }
-        }
         last_print_out = printSteps(last_print_out, start, num_walker, mol_i);
     }
 
-    // write to files
+    experiment.save(motion, dist);
-    const auto folders = make_directory(motion, dist);
-    save_parameter_to_file(std::string("ste"), folders, parameter, optional);
-    save_data_to_file(std::string("coscos"), folders, mixing_times, cc_dict, optional);
-    save_data_to_file(std::string("sinsin"), folders, mixing_times, ss_dict, optional);
-    save_data_to_file(std::string("f2"), folders, mixing_times, f2, optional);
-
     printEnd(start);
 }
 
@@ -174,7 +39,6 @@ Trajectory make_trajectory(
     const double t_max,
     std::mt19937_64& rng
     ) {
-    // Starting position
     double t_passed = 0;
     double phase = 0;
 

src/sims.h

@@ -1,23 +1,22 @@
 #ifndef RWSIM_SIMS_H
 #define RWSIM_SIMS_H
 
+#include "experiments/base.h"
 #include "motions/base.h"
 #include "times/base.h"
 
 #include <unordered_map>
 #include <string>
 #include <chrono>
-#include <vector>
+#include <random>
 
-struct Trajectory {
+void run_simulation(
-    std::vector<double> time;
+    Experiment& experiment,
-    std::vector<double> phase;
+    std::unordered_map<std::string, double>& parameter,
-    std::vector<double> omega;
+    std::unordered_map<std::string, double>& optional,
-};
+    motions::BaseMotion& motion,
-
+    times::BaseDistribution& dist,
-void run_spectrum(std::unordered_map<std::string, double>& parameter, std::unordered_map<std::string, double> optional, motions::BaseMotion& motion, times::BaseDistribution& dist, std::mt19937_64& rng);
+    std::mt19937_64& rng);
-
-void run_ste(std::unordered_map<std::string, double>& parameter, std::unordered_map<std::string, double> optional, motions::BaseMotion& motion, times::BaseDistribution& dist, std::mt19937_64& rng);
 
 Trajectory make_trajectory(motions::BaseMotion& motion, times::BaseDistribution& dist, double t_max, std::mt19937_64& rng);