From e45ef8162d9656c38c1d4df811f566c8ae49ea0b Mon Sep 17 00:00:00 2001 From: Dominik Demuth Date: Sun, 8 Mar 2026 11:43:47 +0100 Subject: [PATCH] Create Experiment class --- src/CMakeLists.txt | 5 +- src/experiments/CMakeLists.txt | 10 +++ src/experiments/base.h | 31 +++++++ src/experiments/spectrum.cpp | 65 ++++++++++++++ src/experiments/spectrum.h | 30 +++++++ src/experiments/ste.cpp | 99 +++++++++++++++++++++ src/experiments/ste.h | 33 +++++++ src/main.cpp | 9 +- src/sims.cpp | 152 ++------------------------------- src/sims.h | 19 ++--- 10 files changed, 295 insertions(+), 158 deletions(-) create mode 100644 src/experiments/CMakeLists.txt create mode 100644 src/experiments/base.h create mode 100644 src/experiments/spectrum.cpp create mode 100644 src/experiments/spectrum.h create mode 100644 src/experiments/ste.cpp create mode 100644 src/experiments/ste.h diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 2f3621e..831339b 100644 --- a/src/CMakeLists.txt +++ b/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) diff --git a/src/experiments/CMakeLists.txt b/src/experiments/CMakeLists.txt new file mode 100644 index 0000000..3958c18 --- /dev/null +++ b/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) diff --git a/src/experiments/base.h b/src/experiments/base.h new file mode 100644 index 0000000..2688545 --- /dev/null +++ b/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 +#include +#include +#include + +struct Trajectory { + std::vector time; + std::vector phase; + std::vector omega; +}; + +class Experiment { +public: + virtual ~Experiment() = default; + + virtual void setup(const std::unordered_map& parameter, + const std::unordered_map& 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 clone() const = 0; + virtual void merge(const Experiment& other) = 0; +}; + +#endif //RWSIM_EXPERIMENTBASE_H diff --git a/src/experiments/spectrum.cpp b/src/experiments/spectrum.cpp new file mode 100644 index 0000000..8477ce6 --- /dev/null +++ b/src/experiments/spectrum.cpp @@ -0,0 +1,65 @@ +#include "spectrum.h" +#include "../utils/functions.h" +#include "../utils/ranges.h" +#include "../utils/io.h" + +#include +#include + +void SpectrumExperiment::setup( + const std::unordered_map& parameter, + const std::unordered_map& optional + ) { + m_parameter = parameter; + m_optional = optional; + m_num_acq = static_cast(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(parameter.at("techo_steps"))); + + m_fid_dict.clear(); + for (auto t_echo_i: m_echo_times) { + m_fid_dict[t_echo_i] = std::vector(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 SpectrumExperiment::clone() const { + return std::make_unique(*this); +} + +void SpectrumExperiment::merge(const Experiment& other) { + const auto& o = dynamic_cast(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]; + } + } +} diff --git a/src/experiments/spectrum.h b/src/experiments/spectrum.h new file mode 100644 index 0000000..2ac40f2 --- /dev/null +++ b/src/experiments/spectrum.h @@ -0,0 +1,30 @@ +#ifndef RWSIM_SPECTRUM_H +#define RWSIM_SPECTRUM_H + +#include "base.h" + +#include +#include +#include + +class SpectrumExperiment final : public Experiment { +public: + void setup(const std::unordered_map& parameter, + const std::unordered_map& 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 clone() const override; + void merge(const Experiment& other) override; + +private: + int m_num_acq{0}; + std::vector m_t_fid; + std::vector m_echo_times; + std::map> m_fid_dict; + std::unordered_map m_parameter; + std::unordered_map m_optional; + double m_tmax{0}; +}; + +#endif //RWSIM_SPECTRUM_H diff --git a/src/experiments/ste.cpp b/src/experiments/ste.cpp new file mode 100644 index 0000000..d1224df --- /dev/null +++ b/src/experiments/ste.cpp @@ -0,0 +1,99 @@ +#include "ste.h" +#include "../utils/functions.h" +#include "../utils/ranges.h" +#include "../utils/io.h" + +#include +#include + +void STEExperiment::setup( + const std::unordered_map& parameter, + const std::unordered_map& optional + ) { + m_parameter = parameter; + m_optional = optional; + m_num_mix_times = static_cast(parameter.at("tmix_steps")); + m_evolution_times = linspace(parameter.at("tevo_start"), parameter.at("tevo_stop"), static_cast(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(m_num_mix_times, 0.); + m_ss_dict[t_evo_i] = std::vector(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 STEExperiment::clone() const { + return std::make_unique(*this); +} + +void STEExperiment::merge(const Experiment& other) { + const auto& o = dynamic_cast(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]; + } +} diff --git a/src/experiments/ste.h b/src/experiments/ste.h new file mode 100644 index 0000000..a276d77 --- /dev/null +++ b/src/experiments/ste.h @@ -0,0 +1,33 @@ +#ifndef RWSIM_STE_H +#define RWSIM_STE_H + +#include "base.h" + +#include +#include +#include + +class STEExperiment final : public Experiment { +public: + void setup(const std::unordered_map& parameter, + const std::unordered_map& 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 clone() const override; + void merge(const Experiment& other) override; + +private: + int m_num_mix_times{0}; + double m_tpulse4{0}; + std::vector m_evolution_times; + std::vector m_mixing_times; + std::map> m_cc_dict; + std::map> m_ss_dict; + std::vector m_f2; + std::unordered_map m_parameter; + std::unordered_map m_optional; + double m_tmax{0}; +}; + +#endif //RWSIM_STE_H diff --git a/src/main.cpp b/src/main.cpp index e415fb0..02861a7 100644 --- a/src/main.cpp +++ b/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; diff --git a/src/sims.cpp b/src/sims.cpp index f07e824..701c285 100644 --- a/src/sims.cpp +++ b/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 -#include -#include -#include -#include -#include -#include #include - -void run_spectrum( +void run_simulation( + Experiment& experiment, std::unordered_map& parameter, - std::unordered_map optional, + std::unordered_map& optional, motions::BaseMotion& motion, times::BaseDistribution& dist, std::mt19937_64& rng ) { const int num_walker = static_cast(parameter["num_walker"]); - // time axis for all time signals - const int num_acq = static_cast(parameter["num_acq"]); - const std::vector t_fid = arange(num_acq, parameter["dwell_time"]); - const std::vector echo_times = linspace(parameter["techo_start"], parameter["techo_stop"], static_cast(parameter["techo_steps"])); - - // make timesignal vectors and set them to zero - std::map> fid_dict; - for (auto t_echo_i: echo_times) { - fid_dict[t_echo_i] = std::vector(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); - - 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& parameter, - std::unordered_map optional, - motions::BaseMotion& motion, - times::BaseDistribution& dist, - std::mt19937_64& rng - ) { - const int num_walker = static_cast(parameter[std::string("num_walker")]); - - const int num_mix_times = static_cast(parameter["tmix_steps"]); - const std::vector evolution_times = linspace(parameter["tevo_start"], parameter["tevo_stop"], static_cast(parameter["tevo_steps"])); - const std::vector 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> cc_dict; - std::map> ss_dict; - for (auto t_evo_i: evolution_times) { - cc_dict[t_evo_i] = std::vector(num_mix_times); - ss_dict[t_evo_i] = std::vector(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 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; - } - } + for (int mol_i = 0; mol_i < num_walker; mol_i++) { + auto traj = make_trajectory(motion, dist, experiment.tmax(), rng); + experiment.accumulate(traj, motion.getInitOmega(), num_walker); last_print_out = printSteps(last_print_out, start, num_walker, mol_i); } - // write to files - 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); - + experiment.save(motion, dist); 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; diff --git a/src/sims.h b/src/sims.h index 1c5fb4b..e64d319 100644 --- a/src/sims.h +++ b/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 #include #include -#include +#include -struct Trajectory { - std::vector time; - std::vector phase; - std::vector omega; -}; - -void run_spectrum(std::unordered_map& parameter, std::unordered_map optional, motions::BaseMotion& motion, times::BaseDistribution& dist, std::mt19937_64& rng); - -void run_ste(std::unordered_map& parameter, std::unordered_map optional, motions::BaseMotion& motion, times::BaseDistribution& dist, std::mt19937_64& rng); +void run_simulation( + Experiment& experiment, + std::unordered_map& parameter, + std::unordered_map& 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);