From 285c78bed54e2144ae7b86b15cb5c10c595a15a1 Mon Sep 17 00:00:00 2001 From: Dominik Demuth Date: Sun, 8 Mar 2026 14:01:37 +0100 Subject: [PATCH] Formatting --- src/CMakeLists.txt | 26 ++- src/dynamics/CMakeLists.txt | 8 +- src/dynamics/base.h | 21 +- src/dynamics/decoupled.cpp | 29 +-- src/dynamics/decoupled.h | 25 +-- src/experiments/CMakeLists.txt | 11 +- src/experiments/base.h | 28 +-- src/experiments/spectrum.cpp | 96 +++++---- src/experiments/spectrum.h | 33 +-- src/experiments/ste.cpp | 164 +++++++------- src/experiments/ste.h | 39 ++-- src/main.cpp | 83 ++++--- src/motions/CMakeLists.txt | 40 ++-- src/motions/base.cpp | 85 ++++---- src/motions/base.h | 63 +++--- src/motions/bimodalangle.cpp | 82 +++---- src/motions/bimodalangle.h | 33 +-- src/motions/conemotion.cpp | 29 +-- src/motions/conemotion.h | 25 +-- src/motions/conewobble.cpp | 37 ++-- src/motions/conewobble.h | 20 +- src/motions/coordinates.cpp | 56 ++--- src/motions/coordinates.h | 28 +-- src/motions/diffusivemotion.cpp | 8 +- src/motions/diffusivemotion.h | 18 +- src/motions/foursitejump.cpp | 62 +++--- src/motions/foursitejump.h | 35 ++- src/motions/isosmallangle.cpp | 62 +++--- src/motions/isosmallangle.h | 29 +-- src/motions/nsiteconejump.cpp | 84 ++++---- src/motions/nsiteconejump.h | 41 ++-- src/motions/random.cpp | 32 ++- src/motions/random.h | 23 +- src/motions/rjoac.cpp | 40 ++-- src/motions/rjoac.h | 20 +- src/motions/sixsitejump.cpp | 68 +++--- src/motions/sixsitejump.h | 36 ++-- src/sims.cpp | 189 ++++++++-------- src/sims.h | 22 +- src/times/CMakeLists.txt | 12 +- src/times/base.cpp | 44 ++-- src/times/base.h | 46 ++-- src/times/delta.cpp | 27 +-- src/times/delta.h | 20 +- src/times/lognormal.cpp | 60 +++--- src/times/lognormal.h | 31 +-- src/utils/functions.cpp | 71 +++--- src/utils/functions.h | 12 +- src/utils/io.cpp | 368 ++++++++++++++++---------------- src/utils/io.h | 39 ++-- src/utils/ranges.cpp | 87 ++++---- src/utils/registry.h | 62 +++--- 52 files changed, 1381 insertions(+), 1328 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 5f861d5..ad3bbbb 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -5,14 +5,24 @@ add_subdirectory(utils) add_subdirectory(experiments) add_subdirectory(dynamics) -add_library(simulation STATIC sims.cpp sims.h) -target_link_libraries(simulation PRIVATE utils experiments dynamics OpenMP::OpenMP_CXX) - -add_executable( - rwsim - main.cpp +add_library(simulation STATIC + sims.cpp sims.h +) +target_link_libraries(simulation PRIVATE + utils + experiments + dynamics + OpenMP::OpenMP_CXX ) -target_link_libraries(rwsim PUBLIC RWMotion RWTime utils experiments dynamics simulation) +add_executable(rwsim main.cpp) + +target_link_libraries(rwsim PUBLIC + RWMotion + RWTime + utils + experiments + dynamics + simulation +) -target_compile_options(rwsim PUBLIC -Werror -Wall -Wextra -Wconversion -O2) diff --git a/src/dynamics/CMakeLists.txt b/src/dynamics/CMakeLists.txt index ec675d2..8a47fcb 100644 --- a/src/dynamics/CMakeLists.txt +++ b/src/dynamics/CMakeLists.txt @@ -1,6 +1,4 @@ -add_library( - dynamics STATIC - base.h - decoupled.cpp - decoupled.h +add_library(dynamics STATIC + base.h + decoupled.cpp decoupled.h ) diff --git a/src/dynamics/base.h b/src/dynamics/base.h index a71e62c..4cf398e 100644 --- a/src/dynamics/base.h +++ b/src/dynamics/base.h @@ -7,20 +7,21 @@ #include struct Step { - double dt; - double omega; + double dt; + double omega; }; class Dynamics { public: - virtual ~Dynamics() = default; + virtual ~Dynamics() = default; - virtual void setParameters(const std::unordered_map& parameters) = 0; - virtual void initialize(std::mt19937_64& rng) = 0; - virtual Step next(std::mt19937_64& rng) = 0; - [[nodiscard]] virtual double getInitOmega() const = 0; - [[nodiscard]] virtual std::unique_ptr clone() const = 0; - [[nodiscard]] virtual std::string toString() const = 0; + virtual void + setParameters(const std::unordered_map ¶meters) = 0; + virtual void initialize(std::mt19937_64 &rng) = 0; + virtual Step next(std::mt19937_64 &rng) = 0; + [[nodiscard]] virtual double getInitOmega() const = 0; + [[nodiscard]] virtual std::unique_ptr clone() const = 0; + [[nodiscard]] virtual std::string toString() const = 0; }; -#endif //RWSIM_DYNAMICSBASE_H +#endif // RWSIM_DYNAMICSBASE_H diff --git a/src/dynamics/decoupled.cpp b/src/dynamics/decoupled.cpp index 681c4cb..d2de017 100644 --- a/src/dynamics/decoupled.cpp +++ b/src/dynamics/decoupled.cpp @@ -5,31 +5,32 @@ DecoupledDynamics::DecoupledDynamics( std::unique_ptr dist) : m_motion(std::move(motion)), m_dist(std::move(dist)) {} -void DecoupledDynamics::setParameters(const std::unordered_map& parameters) { - m_motion->setParameters(parameters); - m_dist->setParameters(parameters); +void DecoupledDynamics::setParameters( + const std::unordered_map ¶meters) { + m_motion->setParameters(parameters); + m_dist->setParameters(parameters); } -void DecoupledDynamics::initialize(std::mt19937_64& rng) { - m_motion->initialize(rng); - m_dist->initialize(rng); +void DecoupledDynamics::initialize(std::mt19937_64 &rng) { + m_motion->initialize(rng); + m_dist->initialize(rng); } -Step DecoupledDynamics::next(std::mt19937_64& rng) { - double dt = m_dist->tau_wait(rng); - double omega = m_motion->jump(rng); - return {dt, omega}; +Step DecoupledDynamics::next(std::mt19937_64 &rng) { + double dt = m_dist->tau_wait(rng); + double omega = m_motion->jump(rng); + return {dt, omega}; } double DecoupledDynamics::getInitOmega() const { - return m_motion->getInitOmega(); + return m_motion->getInitOmega(); } std::unique_ptr DecoupledDynamics::clone() const { - return std::make_unique( - m_motion->clone(), m_dist->clone()); + return std::make_unique(m_motion->clone(), + m_dist->clone()); } std::string DecoupledDynamics::toString() const { - return m_motion->toString() + "/" + m_dist->toString(); + return m_motion->toString() + "/" + m_dist->toString(); } diff --git a/src/dynamics/decoupled.h b/src/dynamics/decoupled.h index 6be181d..81d4d7a 100644 --- a/src/dynamics/decoupled.h +++ b/src/dynamics/decoupled.h @@ -1,27 +1,28 @@ #ifndef RWSIM_DECOUPLED_H #define RWSIM_DECOUPLED_H -#include "base.h" #include "../motions/base.h" #include "../times/base.h" +#include "base.h" #include class DecoupledDynamics final : public Dynamics { public: - DecoupledDynamics(std::unique_ptr motion, - std::unique_ptr dist); + DecoupledDynamics(std::unique_ptr motion, + std::unique_ptr dist); - void setParameters(const std::unordered_map& parameters) override; - void initialize(std::mt19937_64& rng) override; - Step next(std::mt19937_64& rng) override; - [[nodiscard]] double getInitOmega() const override; - [[nodiscard]] std::unique_ptr clone() const override; - [[nodiscard]] std::string toString() const override; + void setParameters( + const std::unordered_map ¶meters) override; + void initialize(std::mt19937_64 &rng) override; + Step next(std::mt19937_64 &rng) override; + [[nodiscard]] double getInitOmega() const override; + [[nodiscard]] std::unique_ptr clone() const override; + [[nodiscard]] std::string toString() const override; private: - std::unique_ptr m_motion; - std::unique_ptr m_dist; + std::unique_ptr m_motion; + std::unique_ptr m_dist; }; -#endif //RWSIM_DECOUPLED_H +#endif // RWSIM_DECOUPLED_H diff --git a/src/experiments/CMakeLists.txt b/src/experiments/CMakeLists.txt index 3958c18..2da8a7b 100644 --- a/src/experiments/CMakeLists.txt +++ b/src/experiments/CMakeLists.txt @@ -1,10 +1,7 @@ -add_library( - experiments STATIC - base.h - spectrum.cpp - spectrum.h - ste.cpp - ste.h +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 index d13c34e..118a734 100644 --- a/src/experiments/base.h +++ b/src/experiments/base.h @@ -2,27 +2,29 @@ #define RWSIM_EXPERIMENTBASE_H #include -#include #include +#include #include struct Trajectory { - std::vector time; - std::vector phase; - std::vector omega; + std::vector time; + std::vector phase; + std::vector omega; }; class Experiment { public: - virtual ~Experiment() = default; + 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 std::string& directory) = 0; - [[nodiscard]] virtual std::unique_ptr clone() const = 0; - virtual void merge(const Experiment& other) = 0; + virtual void + setup(const std::unordered_map ¶meter, + 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 std::string &directory) = 0; + [[nodiscard]] virtual std::unique_ptr clone() const = 0; + virtual void merge(const Experiment &other) = 0; }; -#endif //RWSIM_EXPERIMENTBASE_H +#endif // RWSIM_EXPERIMENTBASE_H diff --git a/src/experiments/spectrum.cpp b/src/experiments/spectrum.cpp index 6439db4..4e37f32 100644 --- a/src/experiments/spectrum.cpp +++ b/src/experiments/spectrum.cpp @@ -1,65 +1,69 @@ #include "spectrum.h" #include "../utils/functions.h" -#include "../utils/ranges.h" #include "../utils/io.h" - +#include "../utils/ranges.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"))); + const std::unordered_map ¶meter, + 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_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; } - - 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 std::string& directory) { - save_parameter_to_file(std::string("timesignal"), directory, m_parameter, m_optional); - save_data_to_file(std::string("timesignal"), directory, m_t_fid, m_fid_dict, m_optional); +void SpectrumExperiment::save(const std::string &directory) { + save_parameter_to_file(std::string("timesignal"), directory, m_parameter, + m_optional); + save_data_to_file(std::string("timesignal"), directory, m_t_fid, m_fid_dict, + m_optional); } std::unique_ptr SpectrumExperiment::clone() const { - return std::make_unique(*this); + 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]; - } +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 index 86155c9..22dfd75 100644 --- a/src/experiments/spectrum.h +++ b/src/experiments/spectrum.h @@ -4,27 +4,28 @@ #include "base.h" #include -#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 std::string& directory) override; - [[nodiscard]] std::unique_ptr clone() const override; - void merge(const Experiment& other) override; + void setup(const std::unordered_map ¶meter, + 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 std::string &directory) 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}; + 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 +#endif // RWSIM_SPECTRUM_H diff --git a/src/experiments/ste.cpp b/src/experiments/ste.cpp index 7808ac5..b85bc5a 100644 --- a/src/experiments/ste.cpp +++ b/src/experiments/ste.cpp @@ -1,98 +1,108 @@ #include "ste.h" #include "../utils/functions.h" -#include "../utils/ranges.h" #include "../utils/io.h" +#include "../utils/ranges.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"); + const std::unordered_map ¶meter, + 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_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; + 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; +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::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 std::string& directory) { - save_parameter_to_file(std::string("ste"), directory, m_parameter, m_optional); - save_data_to_file(std::string("coscos"), directory, m_mixing_times, m_cc_dict, m_optional); - save_data_to_file(std::string("sinsin"), directory, m_mixing_times, m_ss_dict, m_optional); - save_data_to_file(std::string("f2"), directory, m_mixing_times, m_f2, m_optional); +void STEExperiment::save(const std::string &directory) { + save_parameter_to_file(std::string("ste"), directory, m_parameter, + m_optional); + save_data_to_file(std::string("coscos"), directory, m_mixing_times, m_cc_dict, + m_optional); + save_data_to_file(std::string("sinsin"), directory, m_mixing_times, m_ss_dict, + m_optional); + save_data_to_file(std::string("f2"), directory, m_mixing_times, m_f2, + m_optional); } std::unique_ptr STEExperiment::clone() const { - return std::make_unique(*this); + 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]; +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 index 09766cc..a4f831a 100644 --- a/src/experiments/ste.h +++ b/src/experiments/ste.h @@ -4,30 +4,31 @@ #include "base.h" #include -#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 std::string& directory) override; - [[nodiscard]] std::unique_ptr clone() const override; - void merge(const Experiment& other) override; + void setup(const std::unordered_map ¶meter, + 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 std::string &directory) 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}; + 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 +#endif // RWSIM_STE_H diff --git a/src/main.cpp b/src/main.cpp index ad53db3..d16ab1c 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,61 +1,58 @@ -#include "sims.h" #include "dynamics/decoupled.h" #include "experiments/spectrum.h" #include "experiments/ste.h" -#include "utils/io.h" #include "motions/base.h" +#include "sims.h" #include "times/base.h" +#include "utils/io.h" #include -#include #include +#include +int main(const int argc, char *argv[]) { + Arguments args; + try { + args = parse_args(argc, argv); + } catch (std::invalid_argument &error) { + std::cerr << error.what() << std::endl; + return 1; + } + std::unordered_map parameter{read_parameter(args.parameter_file)}; + for (const auto &[key, value] : args.optional) { + parameter[key] = value; + } -int main (const int argc, char *argv[]) { - Arguments args; - try { - args = parse_args(argc, argv); - } catch (std::invalid_argument& error) { - std::cerr << error.what() << std::endl; - return 1; - } + // print parameter of simulation to inform user + std::cout << "Found parameter\n"; + for (const auto &[key, value] : parameter) { + std::cout << key << ": " << std::to_string(value) << "\n"; + } + std::cout << std::endl; - std::unordered_map parameter { read_parameter(args.parameter_file) }; + std::mt19937_64 rng; + if (parameter.count("seed")) { + rng.seed(static_cast(parameter.at("seed"))); + } else { + std::random_device rd; + rng.seed(rd()); + } - for (const auto& [key, value]: args.optional) { - parameter[key] = value; - } + DecoupledDynamics dynamics( + motions::BaseMotion::createFromInput(args.motion_type), + times::BaseDistribution::createFromInput(args.distribution_type)); - // print parameter of simulation to inform user - std::cout << "Found parameter\n"; - for (const auto& [key, value]: parameter) { - std::cout << key << ": " << std::to_string(value) << "\n"; - } - std::cout << std::endl; + if (args.spectrum) { + SpectrumExperiment experiment; + run_simulation(experiment, parameter, args.optional, dynamics, rng); + } + if (args.ste) { + STEExperiment experiment; + run_simulation(experiment, parameter, args.optional, dynamics, rng); + } - std::mt19937_64 rng; - if (parameter.count("seed")) { - rng.seed(static_cast(parameter.at("seed"))); - } else { - std::random_device rd; - rng.seed(rd()); - } - - DecoupledDynamics dynamics( - motions::BaseMotion::createFromInput(args.motion_type), - times::BaseDistribution::createFromInput(args.distribution_type)); - - if (args.spectrum) { - SpectrumExperiment experiment; - run_simulation(experiment, parameter, args.optional, dynamics, rng); - } - if (args.ste) { - STEExperiment experiment; - run_simulation(experiment, parameter, args.optional, dynamics, rng); - } - - return 0; + return 0; } diff --git a/src/motions/CMakeLists.txt b/src/motions/CMakeLists.txt index aa432b8..fd34bce 100644 --- a/src/motions/CMakeLists.txt +++ b/src/motions/CMakeLists.txt @@ -1,29 +1,15 @@ -# Create a library target for motions -add_library( - RWMotion OBJECT - conewobble.cpp - conewobble.h - conemotion.cpp - conemotion.h - diffusivemotion.cpp - diffusivemotion.h - coordinates.cpp - coordinates.h - base.cpp - base.h - random.cpp - random.h - isosmallangle.cpp - isosmallangle.h - foursitejump.cpp - foursitejump.h - rjoac.cpp - rjoac.h - bimodalangle.cpp - bimodalangle.h - sixsitejump.cpp - sixsitejump.h - nsiteconejump.cpp - nsiteconejump.h +add_library(RWMotion OBJECT + base.cpp base.h + bimodalangle.cpp bimodalangle.h + conewobble.cpp conewobble.h + conemotion.cpp conemotion.h + coordinates.cpp coordinates.h + diffusivemotion.cpp diffusivemotion.h + foursitejump.cpp foursitejump.h + isosmallangle.cpp isosmallangle.h + nsiteconejump.cpp nsiteconejump.h + random.cpp random.h + rjoac.cpp rjoac.h + sixsitejump.cpp sixsitejump.h ) diff --git a/src/motions/base.cpp b/src/motions/base.cpp index 7cabf7c..ff81e8a 100644 --- a/src/motions/base.cpp +++ b/src/motions/base.cpp @@ -5,48 +5,51 @@ #include namespace motions { - BaseMotion::BaseMotion(std::string name, const double delta, const double eta) : m_name(std::move(name)), m_delta(delta), m_eta(eta) {} +BaseMotion::BaseMotion(std::string name, const double delta, const double eta) + : m_name(std::move(name)), m_delta(delta), m_eta(eta) {} - BaseMotion::BaseMotion(std::string name) : m_name(std::move(name)) {} +BaseMotion::BaseMotion(std::string name) : m_name(std::move(name)) {} - double BaseMotion::omega_q(const double cos_theta, const double phi) const { - const double cos_theta_square = cos_theta * cos_theta; - const double sin_theta_square = 1. - cos_theta_square; +double BaseMotion::omega_q(const double cos_theta, const double phi) const { + const double cos_theta_square = cos_theta * cos_theta; + const double sin_theta_square = 1. - cos_theta_square; - return M_PI * m_delta * (3. * cos_theta_square - 1. - m_eta * sin_theta_square * std::cos(2.*phi)); - } - - double BaseMotion::omega_q(const coordinates::SphericalPos& pos) const { - auto [cos_theta, phi] = pos; - - return omega_q(cos_theta, phi); - } - - coordinates::SphericalPos BaseMotion::draw_position(std::mt19937_64& rng) { - const double cos_theta = 1 - 2 * m_uni_dist(rng); - const double phi = 2.0 * M_PI * m_uni_dist(rng); - - return {cos_theta, phi}; - } - - std::unique_ptr BaseMotion::createFromInput(const std::string& input) { - return Registry::create(input); - } - - void BaseMotion::setParameters(const std::unordered_map ¶meters) { - m_delta = parameters.at("delta"); - m_eta = parameters.at("eta"); - } - - std::unordered_map BaseMotion::getParameters() const { - return std::unordered_map{ - {"delta", m_delta}, - {"eta", m_eta} - }; - } - - std::ostream& operator<<(std::ostream& os, const BaseMotion& m) { - os << m.getName(); - return os; - } + return M_PI * m_delta * + (3. * cos_theta_square - 1. - + m_eta * sin_theta_square * std::cos(2. * phi)); } + +double BaseMotion::omega_q(const coordinates::SphericalPos &pos) const { + auto [cos_theta, phi] = pos; + + return omega_q(cos_theta, phi); +} + +coordinates::SphericalPos BaseMotion::draw_position(std::mt19937_64 &rng) { + const double cos_theta = 1 - 2 * m_uni_dist(rng); + const double phi = 2.0 * M_PI * m_uni_dist(rng); + + return {cos_theta, phi}; +} + +std::unique_ptr +BaseMotion::createFromInput(const std::string &input) { + return Registry::create(input); +} + +void BaseMotion::setParameters( + const std::unordered_map ¶meters) { + m_delta = parameters.at("delta"); + m_eta = parameters.at("eta"); +} + +std::unordered_map BaseMotion::getParameters() const { + return std::unordered_map{{"delta", m_delta}, + {"eta", m_eta}}; +} + +std::ostream &operator<<(std::ostream &os, const BaseMotion &m) { + os << m.getName(); + return os; +} +} // namespace motions diff --git a/src/motions/base.h b/src/motions/base.h index cc37f23..91500e1 100644 --- a/src/motions/base.h +++ b/src/motions/base.h @@ -8,42 +8,43 @@ #include namespace motions { - class BaseMotion { - public: - virtual ~BaseMotion() = default; +class BaseMotion { +public: + virtual ~BaseMotion() = default; - BaseMotion(std::string, double, double); - explicit BaseMotion(std::string); + BaseMotion(std::string, double, double); + explicit BaseMotion(std::string); - coordinates::SphericalPos draw_position(std::mt19937_64& rng); - [[nodiscard]] double omega_q(double, double) const; - [[nodiscard]] double omega_q(const coordinates::SphericalPos&) const; + coordinates::SphericalPos draw_position(std::mt19937_64 &rng); + [[nodiscard]] double omega_q(double, double) const; + [[nodiscard]] double omega_q(const coordinates::SphericalPos &) const; - virtual void initialize(std::mt19937_64& rng) = 0; - virtual double jump(std::mt19937_64& rng) = 0; - [[nodiscard]] virtual std::unique_ptr clone() const = 0; + virtual void initialize(std::mt19937_64 &rng) = 0; + virtual double jump(std::mt19937_64 &rng) = 0; + [[nodiscard]] virtual std::unique_ptr clone() const = 0; - virtual void setParameters(const std::unordered_map&); - [[nodiscard]] virtual std::unordered_map getParameters() const; - [[nodiscard]] double getDelta() const { return m_delta; } - void setDelta(const double delta) { m_delta = delta; } - [[nodiscard]] double getEta() const { return m_eta; } - void setEta(const double eta) { m_eta = eta; } - [[nodiscard]] std::string getName() const { return m_name; } - [[nodiscard]] double getInitOmega() const { return m_initial_omega; } + virtual void setParameters(const std::unordered_map &); + [[nodiscard]] virtual std::unordered_map + getParameters() const; + [[nodiscard]] double getDelta() const { return m_delta; } + void setDelta(const double delta) { m_delta = delta; } + [[nodiscard]] double getEta() const { return m_eta; } + void setEta(const double eta) { m_eta = eta; } + [[nodiscard]] std::string getName() const { return m_name; } + [[nodiscard]] double getInitOmega() const { return m_initial_omega; } - [[nodiscard]] virtual std::string toString() const = 0; + [[nodiscard]] virtual std::string toString() const = 0; - static std::unique_ptr createFromInput(const std::string& input); + static std::unique_ptr createFromInput(const std::string &input); - protected: - std::string m_name{"BaseMotion"}; - double m_delta{1.}; - double m_eta{0.}; - std::uniform_real_distribution<> m_uni_dist{0., 1.}; - double m_initial_omega{0.}; - }; +protected: + std::string m_name{"BaseMotion"}; + double m_delta{1.}; + double m_eta{0.}; + std::uniform_real_distribution<> m_uni_dist{0., 1.}; + double m_initial_omega{0.}; +}; - std::ostream& operator<<(std::ostream& os, const BaseMotion& m); -} -#endif //RWSIM_MOTIONBASE_H +std::ostream &operator<<(std::ostream &os, const BaseMotion &m); +} // namespace motions +#endif // RWSIM_MOTIONBASE_H diff --git a/src/motions/bimodalangle.cpp b/src/motions/bimodalangle.cpp index 0ecbcae..e480db5 100644 --- a/src/motions/bimodalangle.cpp +++ b/src/motions/bimodalangle.cpp @@ -1,47 +1,53 @@ #include "bimodalangle.h" -#include "coordinates.h" #include "../utils/registry.h" +#include "coordinates.h" -static AutoRegister reg("BimodalAngle"); +static AutoRegister + reg("BimodalAngle"); namespace motions { - BimodalAngle::BimodalAngle(const double delta, const double eta, const double angle1, const double angle2, const double prob) : - DiffusiveMotion(std::string("BimodalAngle"), delta, eta), - m_angle1(angle1 * M_PI / 180.0), - m_angle2(angle2 * M_PI / 180.0), - m_prob(prob) {} - BimodalAngle::BimodalAngle() : DiffusiveMotion(std::string("BimodalAngle")) {} +BimodalAngle::BimodalAngle(const double delta, const double eta, + const double angle1, const double angle2, + const double prob) + : DiffusiveMotion(std::string("BimodalAngle"), delta, eta), + m_angle1(angle1 * M_PI / 180.0), m_angle2(angle2 * M_PI / 180.0), + m_prob(prob) {} +BimodalAngle::BimodalAngle() : DiffusiveMotion(std::string("BimodalAngle")) {} - double BimodalAngle::jump(std::mt19937_64& rng) { - const double angle = m_uni_dist(rng) < m_prob ? m_angle1 : m_angle2; - const double gamma{2 * M_PI * m_uni_dist(rng)}; - m_prev_pos = rotate(m_prev_pos, angle, gamma); +double BimodalAngle::jump(std::mt19937_64 &rng) { + const double angle = m_uni_dist(rng) < m_prob ? m_angle1 : m_angle2; + const double gamma{2 * M_PI * m_uni_dist(rng)}; + m_prev_pos = rotate(m_prev_pos, angle, gamma); - return omega_q(m_prev_pos); - } - - std::unique_ptr BimodalAngle::clone() const { - return std::make_unique(*this); - } - - void BimodalAngle::setParameters(const std::unordered_map ¶meter) { - BaseMotion::setParameters(parameter); - - m_angle1 = parameter.at("angle1") * M_PI / 180.; - m_angle2 = parameter.at("angle2") * M_PI / 180.; - m_prob = parameter.at("probability1"); - } - - std::unordered_map BimodalAngle::getParameters() const { - auto parameter = BaseMotion::getParameters(); - parameter["angle1"] = m_angle1 * 180 / M_PI; - parameter["angle2"] = m_angle2 * 180 / M_PI; - parameter["probability1"] = m_prob; - return parameter; - } - - std::string BimodalAngle::toString() const { - return std::string{"BimodalAngle/angle1=" + std::to_string(m_angle1 * 180 / M_PI) + "/angle2=" + std::to_string(m_angle2 * 180 / M_PI) + "/probability1=" + std::to_string(m_prob)}; - } + return omega_q(m_prev_pos); } + +std::unique_ptr BimodalAngle::clone() const { + return std::make_unique(*this); +} + +void BimodalAngle::setParameters( + const std::unordered_map ¶meter) { + BaseMotion::setParameters(parameter); + + m_angle1 = parameter.at("angle1") * M_PI / 180.; + m_angle2 = parameter.at("angle2") * M_PI / 180.; + m_prob = parameter.at("probability1"); +} + +std::unordered_map BimodalAngle::getParameters() const { + auto parameter = BaseMotion::getParameters(); + parameter["angle1"] = m_angle1 * 180 / M_PI; + parameter["angle2"] = m_angle2 * 180 / M_PI; + parameter["probability1"] = m_prob; + return parameter; +} + +std::string BimodalAngle::toString() const { + return std::string{ + "BimodalAngle/angle1=" + std::to_string(m_angle1 * 180 / M_PI) + + "/angle2=" + std::to_string(m_angle2 * 180 / M_PI) + + "/probability1=" + std::to_string(m_prob)}; +} +} // namespace motions diff --git a/src/motions/bimodalangle.h b/src/motions/bimodalangle.h index d444b1b..f6b6045 100644 --- a/src/motions/bimodalangle.h +++ b/src/motions/bimodalangle.h @@ -5,21 +5,22 @@ #include "diffusivemotion.h" namespace motions { - class BimodalAngle final : public DiffusiveMotion { - public: - BimodalAngle(double, double, double, double, double); - BimodalAngle(); +class BimodalAngle final : public DiffusiveMotion { +public: + BimodalAngle(double, double, double, double, double); + BimodalAngle(); - double jump(std::mt19937_64& rng) override; - [[nodiscard]] std::unique_ptr clone() const override; - void setParameters(const std::unordered_map &) override; - [[nodiscard]] std::unordered_map getParameters() const override; - [[nodiscard]] std::string toString() const override; + double jump(std::mt19937_64 &rng) override; + [[nodiscard]] std::unique_ptr clone() const override; + void setParameters(const std::unordered_map &) override; + [[nodiscard]] std::unordered_map + getParameters() const override; + [[nodiscard]] std::string toString() const override; - protected: - double m_angle1{0}; - double m_angle2{0}; - double m_prob{0}; - }; -} -#endif //BIMODALANGLE_H +protected: + double m_angle1{0}; + double m_angle2{0}; + double m_prob{0}; +}; +} // namespace motions +#endif // BIMODALANGLE_H diff --git a/src/motions/conemotion.cpp b/src/motions/conemotion.cpp index ce06904..dded28d 100644 --- a/src/motions/conemotion.cpp +++ b/src/motions/conemotion.cpp @@ -1,18 +1,19 @@ #include "conemotion.h" namespace motions { - void ConeMotion::initialize(std::mt19937_64& rng) { - m_axis = draw_position(rng); - } - - void ConeMotion::setParameters(const std::unordered_map ¶meters) { - BaseMotion::setParameters(parameters); - m_angle = parameters.at("angle"); - } - - std::unordered_map ConeMotion::getParameters() const { - auto parameter = BaseMotion::getParameters(); - parameter["angle"] = m_angle; - return parameter; - } +void ConeMotion::initialize(std::mt19937_64 &rng) { + m_axis = draw_position(rng); } + +void ConeMotion::setParameters( + const std::unordered_map ¶meters) { + BaseMotion::setParameters(parameters); + m_angle = parameters.at("angle"); +} + +std::unordered_map ConeMotion::getParameters() const { + auto parameter = BaseMotion::getParameters(); + parameter["angle"] = m_angle; + return parameter; +} +} // namespace motions diff --git a/src/motions/conemotion.h b/src/motions/conemotion.h index 9b51cae..76a8e1b 100644 --- a/src/motions/conemotion.h +++ b/src/motions/conemotion.h @@ -5,19 +5,20 @@ #include "coordinates.h" namespace motions { - class ConeMotion : public BaseMotion { - public: - using BaseMotion::BaseMotion; +class ConeMotion : public BaseMotion { +public: + using BaseMotion::BaseMotion; - void initialize(std::mt19937_64& rng) override; + void initialize(std::mt19937_64 &rng) override; - void setParameters(const std::unordered_map &) override; - [[nodiscard]] std::unordered_map getParameters() const override; + void setParameters(const std::unordered_map &) override; + [[nodiscard]] std::unordered_map + getParameters() const override; - protected: - double m_angle{0}; - coordinates::SphericalPos m_axis{1, 0}; - }; -} +protected: + double m_angle{0}; + coordinates::SphericalPos m_axis{1, 0}; +}; +} // namespace motions -#endif //CONEMOTION_H +#endif // CONEMOTION_H diff --git a/src/motions/conewobble.cpp b/src/motions/conewobble.cpp index ea8a74f..13288ff 100644 --- a/src/motions/conewobble.cpp +++ b/src/motions/conewobble.cpp @@ -1,25 +1,28 @@ #include "conewobble.h" -#include "coordinates.h" #include "../utils/registry.h" +#include "coordinates.h" static AutoRegister reg("ConeWobble"); namespace motions { - WobbleCone::WobbleCone(const double delta, const double eta, const double chi) : ConeMotion("Wobble in Cone", delta, eta) { m_angle = chi; } - WobbleCone::WobbleCone() : ConeMotion("Wobble in Cone") {} - - double WobbleCone::jump(std::mt19937_64& rng) { - const double real_angle = m_uni_dist(rng) * m_angle; - const double phi = 2 * M_PI * m_uni_dist(rng); - return omega_q(rotate(m_axis, real_angle, phi)); - } - - std::unique_ptr WobbleCone::clone() const { - return std::make_unique(*this); - } - - std::string WobbleCone::toString() const { - return std::string("ConeWobble/angle=") + std::to_string(m_angle); - } +WobbleCone::WobbleCone(const double delta, const double eta, const double chi) + : ConeMotion("Wobble in Cone", delta, eta) { + m_angle = chi; } +WobbleCone::WobbleCone() : ConeMotion("Wobble in Cone") {} + +double WobbleCone::jump(std::mt19937_64 &rng) { + const double real_angle = m_uni_dist(rng) * m_angle; + const double phi = 2 * M_PI * m_uni_dist(rng); + return omega_q(rotate(m_axis, real_angle, phi)); +} + +std::unique_ptr WobbleCone::clone() const { + return std::make_unique(*this); +} + +std::string WobbleCone::toString() const { + return std::string("ConeWobble/angle=") + std::to_string(m_angle); +} +} // namespace motions diff --git a/src/motions/conewobble.h b/src/motions/conewobble.h index 6de4342..e3d97ad 100644 --- a/src/motions/conewobble.h +++ b/src/motions/conewobble.h @@ -4,14 +4,14 @@ #include "conemotion.h" namespace motions { - class WobbleCone final: public ConeMotion { - public: - WobbleCone(double, double, double); - WobbleCone(); +class WobbleCone final : public ConeMotion { +public: + WobbleCone(double, double, double); + WobbleCone(); - double jump(std::mt19937_64& rng) override; - [[nodiscard]] std::unique_ptr clone() const override; - [[nodiscard]] std::string toString() const override; - }; -} -#endif //CONEWOBBLE_H + double jump(std::mt19937_64 &rng) override; + [[nodiscard]] std::unique_ptr clone() const override; + [[nodiscard]] std::string toString() const override; +}; +} // namespace motions +#endif // CONEWOBBLE_H diff --git a/src/motions/coordinates.cpp b/src/motions/coordinates.cpp index 0141500..bd4dc6c 100644 --- a/src/motions/coordinates.cpp +++ b/src/motions/coordinates.cpp @@ -1,41 +1,43 @@ #include "coordinates.h" #include -#include namespace coordinates { - SphericalPos rotate(const SphericalPos& old_pos, const double alpha, const double beta) { - const double sin_alpha{std::sin(alpha)}; - const double cos_alpha{std::cos(alpha)}; +SphericalPos rotate(const SphericalPos &old_pos, const double alpha, + const double beta) { + const double sin_alpha{std::sin(alpha)}; + const double cos_alpha{std::cos(alpha)}; - const double sin_beta{std::sin(beta)}; - const double cos_beta{std::cos(beta)}; + const double sin_beta{std::sin(beta)}; + const double cos_beta{std::cos(beta)}; - const double cos_theta{old_pos.cos_theta}; + const double cos_theta{old_pos.cos_theta}; - if (std::abs(cos_theta) == 1) { - return xyz_to_spherical(CartesianPos{cos_alpha * cos_beta, cos_alpha * sin_beta, cos_alpha * cos_theta}); - } + if (std::abs(cos_theta) == 1) { + return xyz_to_spherical(CartesianPos{ + cos_alpha * cos_beta, cos_alpha * sin_beta, cos_alpha * cos_theta}); + } - const double norm{std::sqrt(1 - cos_theta * cos_theta) + 1e-15}; + const double norm{std::sqrt(1 - cos_theta * cos_theta) + 1e-15}; - auto [x, y , z] = spherical_to_xyz(old_pos); + auto [x, y, z] = spherical_to_xyz(old_pos); - const auto new_pos = CartesianPos{ - cos_alpha * x + sin_alpha * (-cos_beta * y - sin_beta * x * z) / norm, - cos_alpha * y + sin_alpha * (cos_beta * x - sin_beta * y * z) / norm, - cos_alpha * z + sin_alpha * sin_beta * norm, - }; + const auto new_pos = CartesianPos{ + cos_alpha * x + sin_alpha * (-cos_beta * y - sin_beta * x * z) / norm, + cos_alpha * y + sin_alpha * (cos_beta * x - sin_beta * y * z) / norm, + cos_alpha * z + sin_alpha * sin_beta * norm, + }; - return xyz_to_spherical(new_pos); - } + return xyz_to_spherical(new_pos); +} - CartesianPos spherical_to_xyz(const SphericalPos& pos) { - const double sin_theta = std::sin(std::acos(pos.cos_theta)); - return {sin_theta * std::cos(pos.phi), sin_theta * std::sin(pos.phi), pos.cos_theta}; - } +CartesianPos spherical_to_xyz(const SphericalPos &pos) { + const double sin_theta = std::sin(std::acos(pos.cos_theta)); + return {sin_theta * std::cos(pos.phi), sin_theta * std::sin(pos.phi), + pos.cos_theta}; +} - SphericalPos xyz_to_spherical(const CartesianPos& pos) { - return {pos.z, std::atan2(pos.y, pos.x)}; - } -} \ No newline at end of file +SphericalPos xyz_to_spherical(const CartesianPos &pos) { + return {pos.z, std::atan2(pos.y, pos.x)}; +} +} // namespace coordinates diff --git a/src/motions/coordinates.h b/src/motions/coordinates.h index 3e4fd2f..1b1dc9c 100644 --- a/src/motions/coordinates.h +++ b/src/motions/coordinates.h @@ -3,20 +3,20 @@ #define COORDINATES_H namespace coordinates { - struct SphericalPos { - double cos_theta; - double phi; - }; +struct SphericalPos { + double cos_theta; + double phi; +}; - struct CartesianPos { - double x; - double y; - double z; - }; +struct CartesianPos { + double x; + double y; + double z; +}; - SphericalPos rotate(const SphericalPos&, double, double); - CartesianPos spherical_to_xyz(const SphericalPos&); - SphericalPos xyz_to_spherical(const CartesianPos&); -} +SphericalPos rotate(const SphericalPos &, double, double); +CartesianPos spherical_to_xyz(const SphericalPos &); +SphericalPos xyz_to_spherical(const CartesianPos &); +} // namespace coordinates -#endif //COORDINATES_H +#endif // COORDINATES_H diff --git a/src/motions/diffusivemotion.cpp b/src/motions/diffusivemotion.cpp index 642babe..57e39ff 100644 --- a/src/motions/diffusivemotion.cpp +++ b/src/motions/diffusivemotion.cpp @@ -1,8 +1,8 @@ #include "diffusivemotion.h" namespace motions { - void DiffusiveMotion::initialize(std::mt19937_64& rng) { - m_prev_pos = draw_position(rng); - m_initial_omega = omega_q(m_prev_pos); - } +void DiffusiveMotion::initialize(std::mt19937_64 &rng) { + m_prev_pos = draw_position(rng); + m_initial_omega = omega_q(m_prev_pos); } +} // namespace motions diff --git a/src/motions/diffusivemotion.h b/src/motions/diffusivemotion.h index a6e9a28..0ec37c2 100644 --- a/src/motions/diffusivemotion.h +++ b/src/motions/diffusivemotion.h @@ -5,15 +5,15 @@ #include "coordinates.h" namespace motions { - class DiffusiveMotion : public BaseMotion { - public: - using BaseMotion::BaseMotion; +class DiffusiveMotion : public BaseMotion { +public: + using BaseMotion::BaseMotion; - void initialize(std::mt19937_64& rng) override; + void initialize(std::mt19937_64 &rng) override; - protected: - coordinates::SphericalPos m_prev_pos{0., 0.}; - }; -} +protected: + coordinates::SphericalPos m_prev_pos{0., 0.}; +}; +} // namespace motions -#endif //DIFFUSIVEMOTION_H +#endif // DIFFUSIVEMOTION_H diff --git a/src/motions/foursitejump.cpp b/src/motions/foursitejump.cpp index b9cdb7e..94ac805 100644 --- a/src/motions/foursitejump.cpp +++ b/src/motions/foursitejump.cpp @@ -1,40 +1,42 @@ #include "foursitejump.h" -#include "coordinates.h" #include "../utils/registry.h" +#include "coordinates.h" -static AutoRegister reg("FourSiteTetrahedral"); +static AutoRegister + reg("FourSiteTetrahedral"); namespace motions { - FourSiteTetrahedron::FourSiteTetrahedron(const double delta, const double eta) : - BaseMotion(std::string{"FourSiteTetrahedral"}, delta, eta) {} +FourSiteTetrahedron::FourSiteTetrahedron(const double delta, const double eta) + : BaseMotion(std::string{"FourSiteTetrahedral"}, delta, eta) {} - FourSiteTetrahedron::FourSiteTetrahedron() : BaseMotion(std::string{"FourSiteTetrahedral"}) {} +FourSiteTetrahedron::FourSiteTetrahedron() + : BaseMotion(std::string{"FourSiteTetrahedral"}) {} - void FourSiteTetrahedron::initialize(std::mt19937_64& rng) { - const auto pos = draw_position(rng); - m_corners[0] = omega_q(pos); +void FourSiteTetrahedron::initialize(std::mt19937_64 &rng) { + const auto pos = draw_position(rng); + m_corners[0] = omega_q(pos); - const double alpha = 2. * M_PI * m_uni_dist(rng); + const double alpha = 2. * M_PI * m_uni_dist(rng); - for (int i = 1; i<4; i++) { - auto corner_pos = rotate(pos, m_beta, alpha + (i-1) * 2*M_PI/3.); - m_corners[i] = omega_q(corner_pos); - } - m_initial_omega = FourSiteTetrahedron::jump(rng); - } - - double FourSiteTetrahedron::jump(std::mt19937_64& rng) { - m_corner_idx += m_chooser(rng); - m_corner_idx %= 4; - - return m_corners[m_corner_idx]; - } - - std::unique_ptr FourSiteTetrahedron::clone() const { - return std::make_unique(*this); - } - - std::string FourSiteTetrahedron::toString() const { - return {"FourSiteTetrahedral"}; - } + for (int i = 1; i < 4; i++) { + auto corner_pos = rotate(pos, m_beta, alpha + (i - 1) * 2 * M_PI / 3.); + m_corners[i] = omega_q(corner_pos); + } + m_initial_omega = FourSiteTetrahedron::jump(rng); } + +double FourSiteTetrahedron::jump(std::mt19937_64 &rng) { + m_corner_idx += m_chooser(rng); + m_corner_idx %= 4; + + return m_corners[m_corner_idx]; +} + +std::unique_ptr FourSiteTetrahedron::clone() const { + return std::make_unique(*this); +} + +std::string FourSiteTetrahedron::toString() const { + return {"FourSiteTetrahedral"}; +} +} // namespace motions diff --git a/src/motions/foursitejump.h b/src/motions/foursitejump.h index 19bc0db..dd31c94 100644 --- a/src/motions/foursitejump.h +++ b/src/motions/foursitejump.h @@ -3,29 +3,28 @@ #include "base.h" -#include #include +#include namespace motions { - class FourSiteTetrahedron final : public BaseMotion { - public: - FourSiteTetrahedron(double, double); - FourSiteTetrahedron(); +class FourSiteTetrahedron final : public BaseMotion { +public: + FourSiteTetrahedron(double, double); + FourSiteTetrahedron(); - void initialize(std::mt19937_64& rng) override; - double jump(std::mt19937_64& rng) override; - [[nodiscard]] std::unique_ptr clone() const override; + void initialize(std::mt19937_64 &rng) override; + double jump(std::mt19937_64 &rng) override; + [[nodiscard]] std::unique_ptr clone() const override; - [[nodiscard]] std::string toString() const override; + [[nodiscard]] std::string toString() const override; - private: - const double m_beta{std::acos(-1/3.)}; +private: + const double m_beta{std::acos(-1 / 3.)}; - std::array m_corners{}; - int m_corner_idx{0}; + std::array m_corners{}; + int m_corner_idx{0}; - std::uniform_int_distribution<> m_chooser{1, 3}; - - }; -} -#endif //RWSIM_MOTIONTETRAHEDRAL_H + std::uniform_int_distribution<> m_chooser{1, 3}; +}; +} // namespace motions +#endif // RWSIM_MOTIONTETRAHEDRAL_H diff --git a/src/motions/isosmallangle.cpp b/src/motions/isosmallangle.cpp index 1afe7e4..8fa3d9c 100644 --- a/src/motions/isosmallangle.cpp +++ b/src/motions/isosmallangle.cpp @@ -1,37 +1,41 @@ #include "isosmallangle.h" -#include "coordinates.h" #include "../utils/registry.h" +#include "coordinates.h" -static AutoRegister reg("IsotropicAngle"); +static AutoRegister + reg("IsotropicAngle"); namespace motions { - SmallAngle::SmallAngle(const double delta, const double eta, const double chi) : - DiffusiveMotion(std::string("IsotropicAngle"), delta, eta), m_chi(chi * M_PI / 180.0) {} - SmallAngle::SmallAngle() : DiffusiveMotion(std::string("IsotropicAngle")) {} +SmallAngle::SmallAngle(const double delta, const double eta, const double chi) + : DiffusiveMotion(std::string("IsotropicAngle"), delta, eta), + m_chi(chi * M_PI / 180.0) {} +SmallAngle::SmallAngle() : DiffusiveMotion(std::string("IsotropicAngle")) {} - double SmallAngle::jump(std::mt19937_64& rng) { - const double gamma{2 * M_PI * m_uni_dist(rng)}; - m_prev_pos = rotate(m_prev_pos, m_chi, gamma); +double SmallAngle::jump(std::mt19937_64 &rng) { + const double gamma{2 * M_PI * m_uni_dist(rng)}; + m_prev_pos = rotate(m_prev_pos, m_chi, gamma); - return omega_q(m_prev_pos); - } - - std::unique_ptr SmallAngle::clone() const { - return std::make_unique(*this); - } - - void SmallAngle::setParameters(const std::unordered_map ¶meters) { - m_chi = parameters.at("angle") * M_PI / 180.0; - BaseMotion::setParameters(parameters); - } - - std::unordered_map SmallAngle::getParameters() const { - auto parameter = BaseMotion::getParameters(); - parameter["angle"] = m_chi * 180 / M_PI; - return parameter; - } - - std::string SmallAngle::toString() const { - return std::string{"IsotropicAngle/angle=" + std::to_string(m_chi * 180 / M_PI)}; - } + return omega_q(m_prev_pos); } + +std::unique_ptr SmallAngle::clone() const { + return std::make_unique(*this); +} + +void SmallAngle::setParameters( + const std::unordered_map ¶meters) { + m_chi = parameters.at("angle") * M_PI / 180.0; + BaseMotion::setParameters(parameters); +} + +std::unordered_map SmallAngle::getParameters() const { + auto parameter = BaseMotion::getParameters(); + parameter["angle"] = m_chi * 180 / M_PI; + return parameter; +} + +std::string SmallAngle::toString() const { + return std::string{"IsotropicAngle/angle=" + + std::to_string(m_chi * 180 / M_PI)}; +} +} // namespace motions diff --git a/src/motions/isosmallangle.h b/src/motions/isosmallangle.h index 6f431c3..e02323e 100644 --- a/src/motions/isosmallangle.h +++ b/src/motions/isosmallangle.h @@ -5,19 +5,20 @@ #include "diffusivemotion.h" namespace motions { - class SmallAngle final : public DiffusiveMotion { - public: - SmallAngle(double, double, double); - SmallAngle(); +class SmallAngle final : public DiffusiveMotion { +public: + SmallAngle(double, double, double); + SmallAngle(); - double jump(std::mt19937_64& rng) override; - [[nodiscard]] std::unique_ptr clone() const override; - void setParameters(const std::unordered_map &) override; - [[nodiscard]] std::unordered_map getParameters() const override; - [[nodiscard]] std::string toString() const override; + double jump(std::mt19937_64 &rng) override; + [[nodiscard]] std::unique_ptr clone() const override; + void setParameters(const std::unordered_map &) override; + [[nodiscard]] std::unordered_map + getParameters() const override; + [[nodiscard]] std::string toString() const override; - private: - double m_chi{0}; - }; -} -#endif //RWSIM_MOTIONISOSMALLANGLE_H +private: + double m_chi{0}; +}; +} // namespace motions +#endif // RWSIM_MOTIONISOSMALLANGLE_H diff --git a/src/motions/nsiteconejump.cpp b/src/motions/nsiteconejump.cpp index cdf88dd..259addb 100644 --- a/src/motions/nsiteconejump.cpp +++ b/src/motions/nsiteconejump.cpp @@ -1,61 +1,63 @@ - #include "nsiteconejump.h" -#include "coordinates.h" #include "../utils/registry.h" +#include "coordinates.h" #include #include -static AutoRegister reg("NSiteConeJump"); +static AutoRegister + reg("NSiteConeJump"); namespace motions { - NSiteJumpOnCone::NSiteJumpOnCone(const double delta, const double eta, const int num_sites, const double chi) : - BaseMotion("NSiteJumpOnCone", delta, eta), - m_num_sites(num_sites), - m_chi(chi*M_PI/180.) {} +NSiteJumpOnCone::NSiteJumpOnCone(const double delta, const double eta, + const int num_sites, const double chi) + : BaseMotion("NSiteJumpOnCone", delta, eta), m_num_sites(num_sites), + m_chi(chi * M_PI / 180.) {} - NSiteJumpOnCone::NSiteJumpOnCone() : BaseMotion("NSiteJumpOnCone") { } +NSiteJumpOnCone::NSiteJumpOnCone() : BaseMotion("NSiteJumpOnCone") {} - void NSiteJumpOnCone::initialize(std::mt19937_64& rng) { - m_sites = std::vector(m_num_sites); - m_chooser = std::uniform_int_distribution<>{1, m_num_sites - 1}; +void NSiteJumpOnCone::initialize(std::mt19937_64 &rng) { + m_sites = std::vector(m_num_sites); + m_chooser = std::uniform_int_distribution<>{1, m_num_sites - 1}; - m_axis = draw_position(rng); + m_axis = draw_position(rng); - const double alpha = m_uni_dist(rng) * 2 * M_PI; - const double steps = 2*M_PI / m_num_sites; - for (int i = 0; i < m_num_sites; i++) { - m_sites[i] = omega_q(rotate(m_axis, m_chi, i * steps + alpha)); - } - } + const double alpha = m_uni_dist(rng) * 2 * M_PI; + const double steps = 2 * M_PI / m_num_sites; + for (int i = 0; i < m_num_sites; i++) { + m_sites[i] = omega_q(rotate(m_axis, m_chi, i * steps + alpha)); + } +} - double NSiteJumpOnCone::jump(std::mt19937_64& rng) { - m_cone_idx += m_chooser(rng); - m_cone_idx %= m_num_sites; +double NSiteJumpOnCone::jump(std::mt19937_64 &rng) { + m_cone_idx += m_chooser(rng); + m_cone_idx %= m_num_sites; - return m_sites[m_cone_idx]; - } + return m_sites[m_cone_idx]; +} - std::unique_ptr NSiteJumpOnCone::clone() const { - return std::make_unique(*this); - } +std::unique_ptr NSiteJumpOnCone::clone() const { + return std::make_unique(*this); +} - void NSiteJumpOnCone::setParameters(const std::unordered_map ¶meters) { - BaseMotion::setParameters(parameters); - m_num_sites = static_cast(parameters.at("num_sites")); - m_chi = parameters.at("chi") * M_PI/180.; - } +void NSiteJumpOnCone::setParameters( + const std::unordered_map ¶meters) { + BaseMotion::setParameters(parameters); + m_num_sites = static_cast(parameters.at("num_sites")); + m_chi = parameters.at("chi") * M_PI / 180.; +} - std::unordered_map NSiteJumpOnCone::getParameters() const { - auto parameter = BaseMotion::getParameters(); - parameter["num_sites"] = m_num_sites; - parameter["chi"] = m_chi * 180. / M_PI; +std::unordered_map NSiteJumpOnCone::getParameters() const { + auto parameter = BaseMotion::getParameters(); + parameter["num_sites"] = m_num_sites; + parameter["chi"] = m_chi * 180. / M_PI; - return parameter; - } + return parameter; +} - std::string NSiteJumpOnCone::toString() const { - return std::to_string(m_num_sites) + "SiteJumpOnCone/angle=" + std::to_string(m_chi*180/M_PI); - } +std::string NSiteJumpOnCone::toString() const { + return std::to_string(m_num_sites) + + "SiteJumpOnCone/angle=" + std::to_string(m_chi * 180 / M_PI); +} -} // motions +} // namespace motions diff --git a/src/motions/nsiteconejump.h b/src/motions/nsiteconejump.h index 3d0f3d7..a99cae2 100644 --- a/src/motions/nsiteconejump.h +++ b/src/motions/nsiteconejump.h @@ -6,28 +6,29 @@ #include namespace motions { - class NSiteJumpOnCone final : public BaseMotion { - public: - NSiteJumpOnCone(double, double, int, double); - NSiteJumpOnCone(); +class NSiteJumpOnCone final : public BaseMotion { +public: + NSiteJumpOnCone(double, double, int, double); + NSiteJumpOnCone(); - void initialize(std::mt19937_64& rng) override; - double jump(std::mt19937_64& rng) override; - [[nodiscard]] std::unique_ptr clone() const override; + void initialize(std::mt19937_64 &rng) override; + double jump(std::mt19937_64 &rng) override; + [[nodiscard]] std::unique_ptr clone() const override; - [[nodiscard]] std::string toString() const override; + [[nodiscard]] std::string toString() const override; - void setParameters(const std::unordered_map &) override; - [[nodiscard]] std::unordered_map getParameters() const override; + void setParameters(const std::unordered_map &) override; + [[nodiscard]] std::unordered_map + getParameters() const override; - private: - int m_num_sites{2}; - std::vector m_sites{}; - double m_chi{1./2.}; - int m_cone_idx = 0; - coordinates::SphericalPos m_axis{1, 0}; - std::uniform_int_distribution<> m_chooser; - }; -} // motions +private: + int m_num_sites{2}; + std::vector m_sites{}; + double m_chi{1. / 2.}; + int m_cone_idx = 0; + coordinates::SphericalPos m_axis{1, 0}; + std::uniform_int_distribution<> m_chooser; +}; +} // namespace motions -#endif //NSITEJUMPONCONE_H +#endif // NSITEJUMPONCONE_H diff --git a/src/motions/random.cpp b/src/motions/random.cpp index 5719a85..d4c8348 100644 --- a/src/motions/random.cpp +++ b/src/motions/random.cpp @@ -1,27 +1,25 @@ - #include "random.h" #include "../utils/registry.h" static AutoRegister reg("RandomJump"); namespace motions { - RandomJump::RandomJump(const double delta, const double eta) : BaseMotion(std::string("RandomJump"), delta, eta) {} +RandomJump::RandomJump(const double delta, const double eta) + : BaseMotion(std::string("RandomJump"), delta, eta) {} - RandomJump::RandomJump() : BaseMotion(std::string("RandomJump")) {} +RandomJump::RandomJump() : BaseMotion(std::string("RandomJump")) {} - std::string RandomJump::toString() const { - return {"RandomJump"}; - } +std::string RandomJump::toString() const { return {"RandomJump"}; } - std::unique_ptr RandomJump::clone() const { - return std::make_unique(*this); - } - - void RandomJump::initialize(std::mt19937_64& rng) { - m_initial_omega = RandomJump::jump(rng); - } - - double RandomJump::jump(std::mt19937_64& rng) { - return omega_q(draw_position(rng)); - } +std::unique_ptr RandomJump::clone() const { + return std::make_unique(*this); } + +void RandomJump::initialize(std::mt19937_64 &rng) { + m_initial_omega = RandomJump::jump(rng); +} + +double RandomJump::jump(std::mt19937_64 &rng) { + return omega_q(draw_position(rng)); +} +} // namespace motions diff --git a/src/motions/random.h b/src/motions/random.h index ea86246..3a0ef55 100644 --- a/src/motions/random.h +++ b/src/motions/random.h @@ -1,21 +1,20 @@ #ifndef RWSIM_MOTIONRANDOMJUMP_H #define RWSIM_MOTIONRANDOMJUMP_H - #include "base.h" namespace motions { - class RandomJump final : public BaseMotion { - public: - RandomJump(double, double); - RandomJump(); +class RandomJump final : public BaseMotion { +public: + RandomJump(double, double); + RandomJump(); - [[nodiscard]] std::string toString() const override; - [[nodiscard]] std::unique_ptr clone() const override; + [[nodiscard]] std::string toString() const override; + [[nodiscard]] std::unique_ptr clone() const override; - void initialize(std::mt19937_64& rng) override; - double jump(std::mt19937_64& rng) override; - }; -} + void initialize(std::mt19937_64 &rng) override; + double jump(std::mt19937_64 &rng) override; +}; +} // namespace motions -#endif //RWSIM_MOTIONRANDOMJUMP_H +#endif // RWSIM_MOTIONRANDOMJUMP_H diff --git a/src/motions/rjoac.cpp b/src/motions/rjoac.cpp index 6250f68..20aca96 100644 --- a/src/motions/rjoac.cpp +++ b/src/motions/rjoac.cpp @@ -1,24 +1,28 @@ - #include "rjoac.h" -#include "coordinates.h" #include "../utils/registry.h" +#include "coordinates.h" -static AutoRegister reg("RandomJumpOnCone"); +static AutoRegister + reg("RandomJumpOnCone"); namespace motions { - RandomJumpOnCone::RandomJumpOnCone(const double delta, const double eta, const double chi) : ConeMotion("RJ on a Cone", delta, eta) { m_angle = chi; } - RandomJumpOnCone::RandomJumpOnCone() : ConeMotion("RJ on a Cone") {} - - double RandomJumpOnCone::jump(std::mt19937_64& rng) { - const double phi = 2 * M_PI * m_uni_dist(rng); - return omega_q(rotate(m_axis, m_angle, phi)); - } - - std::unique_ptr RandomJumpOnCone::clone() const { - return std::make_unique(*this); - } - - std::string RandomJumpOnCone::toString() const { - return std::string("RandomJumpOnCone/angle=") + std::to_string(m_angle); - } +RandomJumpOnCone::RandomJumpOnCone(const double delta, const double eta, + const double chi) + : ConeMotion("RJ on a Cone", delta, eta) { + m_angle = chi; } +RandomJumpOnCone::RandomJumpOnCone() : ConeMotion("RJ on a Cone") {} + +double RandomJumpOnCone::jump(std::mt19937_64 &rng) { + const double phi = 2 * M_PI * m_uni_dist(rng); + return omega_q(rotate(m_axis, m_angle, phi)); +} + +std::unique_ptr RandomJumpOnCone::clone() const { + return std::make_unique(*this); +} + +std::string RandomJumpOnCone::toString() const { + return std::string("RandomJumpOnCone/angle=") + std::to_string(m_angle); +} +} // namespace motions diff --git a/src/motions/rjoac.h b/src/motions/rjoac.h index fd74ec7..6411a3f 100644 --- a/src/motions/rjoac.h +++ b/src/motions/rjoac.h @@ -4,15 +4,15 @@ #include "conemotion.h" namespace motions { - class RandomJumpOnCone final: public ConeMotion { - public: - RandomJumpOnCone(double, double, double); - RandomJumpOnCone(); +class RandomJumpOnCone final : public ConeMotion { +public: + RandomJumpOnCone(double, double, double); + RandomJumpOnCone(); - double jump(std::mt19937_64& rng) override; - [[nodiscard]] std::unique_ptr clone() const override; - [[nodiscard]] std::string toString() const override; - }; -} + double jump(std::mt19937_64 &rng) override; + [[nodiscard]] std::unique_ptr clone() const override; + [[nodiscard]] std::string toString() const override; +}; +} // namespace motions -#endif //RJOAC_H +#endif // RJOAC_H diff --git a/src/motions/sixsitejump.cpp b/src/motions/sixsitejump.cpp index cafd16e..2675cb2 100644 --- a/src/motions/sixsitejump.cpp +++ b/src/motions/sixsitejump.cpp @@ -1,43 +1,47 @@ #include "sixsitejump.h" -#include "coordinates.h" #include "../utils/registry.h" +#include "coordinates.h" -static AutoRegister reg("SixSiteOctahedralC3"); +static AutoRegister + reg("SixSiteOctahedralC3"); namespace motions { - SixSiteOctahedronC3::SixSiteOctahedronC3(const double delta, const double eta, const double chi) : - BaseMotion(std::string{"SixSiteOctahedral"}, delta, eta), - m_chi{chi*M_PI/180.} {} +SixSiteOctahedronC3::SixSiteOctahedronC3(const double delta, const double eta, + const double chi) + : BaseMotion(std::string{"SixSiteOctahedral"}, delta, eta), + m_chi{chi * M_PI / 180.} {} - SixSiteOctahedronC3::SixSiteOctahedronC3() : BaseMotion(std::string{"SixSiteOctahedralC3"}) {} +SixSiteOctahedronC3::SixSiteOctahedronC3() + : BaseMotion(std::string{"SixSiteOctahedralC3"}) {} - void SixSiteOctahedronC3::initialize(std::mt19937_64& rng) { - const coordinates::SphericalPos c3_axis = draw_position(rng); - const double alpha = 2. * M_PI * m_uni_dist(rng); +void SixSiteOctahedronC3::initialize(std::mt19937_64 &rng) { + const coordinates::SphericalPos c3_axis = draw_position(rng); + const double alpha = 2. * M_PI * m_uni_dist(rng); - const double m_chi_opposite = M_PI - m_chi; + const double m_chi_opposite = M_PI - m_chi; - for (int i = 0; i<3; i++) { - m_corners[2*i] = omega_q(rotate(c3_axis, m_chi, alpha + i * 2./3.*M_PI)); - m_corners[2*i+1] = omega_q(rotate(c3_axis, m_chi_opposite, alpha + i * 2./3.*M_PI + M_PI/3.)); - } + for (int i = 0; i < 3; i++) { + m_corners[2 * i] = + omega_q(rotate(c3_axis, m_chi, alpha + i * 2. / 3. * M_PI)); + m_corners[2 * i + 1] = omega_q(rotate( + c3_axis, m_chi_opposite, alpha + i * 2. / 3. * M_PI + M_PI / 3.)); + } - m_initial_omega = SixSiteOctahedronC3::jump(rng); - } - - - double SixSiteOctahedronC3::jump(std::mt19937_64& rng) { - m_corner_idx += m_chooser(rng); - m_corner_idx %= 6; - - return m_corners[m_corner_idx]; - } - - std::unique_ptr SixSiteOctahedronC3::clone() const { - return std::make_unique(*this); - } - - std::string SixSiteOctahedronC3::toString() const { - return {"SixSiteOctahedral/angle=" + std::to_string(m_chi / M_PI * 180.)}; - } + m_initial_omega = SixSiteOctahedronC3::jump(rng); } + +double SixSiteOctahedronC3::jump(std::mt19937_64 &rng) { + m_corner_idx += m_chooser(rng); + m_corner_idx %= 6; + + return m_corners[m_corner_idx]; +} + +std::unique_ptr SixSiteOctahedronC3::clone() const { + return std::make_unique(*this); +} + +std::string SixSiteOctahedronC3::toString() const { + return {"SixSiteOctahedral/angle=" + std::to_string(m_chi / M_PI * 180.)}; +} +} // namespace motions diff --git a/src/motions/sixsitejump.h b/src/motions/sixsitejump.h index d900b62..b0db52d 100644 --- a/src/motions/sixsitejump.h +++ b/src/motions/sixsitejump.h @@ -1,32 +1,30 @@ - #ifndef SIXSITEJUMP_H #define SIXSITEJUMP_H #include "base.h" -#include #include +#include namespace motions { - class SixSiteOctahedronC3 final : public BaseMotion { - public: - SixSiteOctahedronC3(double, double, double); - SixSiteOctahedronC3(); +class SixSiteOctahedronC3 final : public BaseMotion { +public: + SixSiteOctahedronC3(double, double, double); + SixSiteOctahedronC3(); - void initialize(std::mt19937_64& rng) override; - double jump(std::mt19937_64& rng) override; - [[nodiscard]] std::unique_ptr clone() const override; + void initialize(std::mt19937_64 &rng) override; + double jump(std::mt19937_64 &rng) override; + [[nodiscard]] std::unique_ptr clone() const override; - [[nodiscard]] std::string toString() const override; + [[nodiscard]] std::string toString() const override; - private: - double m_chi{std::acos(-1.0 / 3.0)}; // 54.74 deg +private: + double m_chi{std::acos(-1.0 / 3.0)}; // 54.74 deg - std::array m_corners{}; - int m_corner_idx{0}; + std::array m_corners{}; + int m_corner_idx{0}; - std::uniform_int_distribution<> m_chooser{1, 5}; - - }; -} -#endif //SIXSITEJUMP_H + std::uniform_int_distribution<> m_chooser{1, 5}; +}; +} // namespace motions +#endif // SIXSITEJUMP_H diff --git a/src/sims.cpp b/src/sims.cpp index 97bd071..1c10ec6 100644 --- a/src/sims.cpp +++ b/src/sims.cpp @@ -2,131 +2,124 @@ #include "utils/functions.h" #include "utils/io.h" -#include #include +#include #include +void run_simulation(Experiment &experiment, + const std::unordered_map ¶meter, + const std::unordered_map &optional, + Dynamics &dynamics, std::mt19937_64 &rng) { + const int num_walker = static_cast(parameter.at("num_walker")); -void run_simulation( - Experiment& experiment, - const std::unordered_map& parameter, - const std::unordered_map& optional, - Dynamics& dynamics, - std::mt19937_64& rng - ) { - const int num_walker = static_cast(parameter.at("num_walker")); + dynamics.setParameters(parameter); + experiment.setup(parameter, optional); - dynamics.setParameters(parameter); - experiment.setup(parameter, optional); + const auto start = printStart(optional); - const auto start = printStart(optional); + const int num_threads = omp_get_max_threads(); - const int num_threads = omp_get_max_threads(); + // Create per-thread RNGs seeded deterministically from the main RNG + std::vector thread_rngs; + thread_rngs.reserve(num_threads); + for (int i = 0; i < num_threads; i++) { + thread_rngs.emplace_back(rng()); + } - // Create per-thread RNGs seeded deterministically from the main RNG - std::vector thread_rngs; - thread_rngs.reserve(num_threads); - for (int i = 0; i < num_threads; i++) { - thread_rngs.emplace_back(rng()); + // Create per-thread clones of dynamics and experiment + std::vector> thread_dynamics; + std::vector> thread_experiments; + for (int i = 0; i < num_threads; i++) { + thread_dynamics.push_back(dynamics.clone()); + thread_experiments.push_back(experiment.clone()); + } + + int steps_done = 0; + auto last_print_out = std::chrono::system_clock::now(); + +#pragma omp parallel + { + 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]; + +#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); + + if (tid == 0) { +#pragma omp atomic + steps_done++; + last_print_out = + printSteps(last_print_out, start, num_walker, steps_done); + } else { +#pragma omp atomic + steps_done++; + } } + } - // Create per-thread clones of dynamics and experiment - std::vector> thread_dynamics; - std::vector> thread_experiments; - for (int i = 0; i < num_threads; i++) { - thread_dynamics.push_back(dynamics.clone()); - thread_experiments.push_back(experiment.clone()); - } + // Merge per-thread results + for (int i = 0; i < num_threads; i++) { + experiment.merge(*thread_experiments[i]); + } - int steps_done = 0; - auto last_print_out = std::chrono::system_clock::now(); - - #pragma omp parallel - { - 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]; - - #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); - - if (tid == 0) { - #pragma omp atomic - steps_done++; - last_print_out = printSteps(last_print_out, start, num_walker, steps_done); - } else { - #pragma omp atomic - steps_done++; - } - } - } - - // Merge per-thread results - for (int i = 0; i < num_threads; i++) { - experiment.merge(*thread_experiments[i]); - } - - const auto directory = make_directory(dynamics.toString()); - experiment.save(directory); - printEnd(start); + const auto directory = make_directory(dynamics.toString()); + experiment.save(directory); + printEnd(start); } +Trajectory make_trajectory(Dynamics &dynamics, const double t_max, + std::mt19937_64 &rng) { + double t_passed = 0; + double phase = 0; -Trajectory make_trajectory( - Dynamics& dynamics, - const double t_max, - std::mt19937_64& rng - ) { - double t_passed = 0; - double phase = 0; + dynamics.initialize(rng); - dynamics.initialize(rng); + double omega = dynamics.getInitOmega(); - double omega = dynamics.getInitOmega(); + Trajectory traj; + traj.time.emplace_back(t_passed); + traj.phase.emplace_back(phase); + traj.omega.emplace_back(omega); + + while (t_passed < t_max) { + auto [dt, new_omega] = dynamics.next(rng); + phase += omega * dt; + t_passed += dt; + omega = new_omega; - Trajectory traj; traj.time.emplace_back(t_passed); traj.phase.emplace_back(phase); traj.omega.emplace_back(omega); + } - while (t_passed < t_max) { - auto [dt, new_omega] = dynamics.next(rng); - phase += omega * dt; - t_passed += dt; - omega = new_omega; - - traj.time.emplace_back(t_passed); - traj.phase.emplace_back(phase); - traj.omega.emplace_back(omega); - } - - return traj; + return traj; } +std::chrono::system_clock::time_point +printStart(const std::unordered_map &optional) { + const auto start = std::chrono::system_clock::now(); + const time_t start_time = std::chrono::system_clock::to_time_t(start); -std::chrono::system_clock::time_point printStart(const std::unordered_map &optional) { - const auto start = std::chrono::system_clock::now(); - const time_t start_time = std::chrono::system_clock::to_time_t(start); + std::cout << "Random walk for "; + for (const auto &[key, value] : optional) { + std::cout << key << " = " << value << "; "; + } + std::cout << std::endl; + std::cout << "Start: " << ctime(&start_time); - std::cout << "Random walk for "; - for (const auto& [key, value] : optional) { - std::cout << key << " = " << value << "; "; - } - std::cout << std::endl; - std::cout << "Start: " << ctime(&start_time); - - return start; + return start; } - void printEnd(const std::chrono::system_clock::time_point start) { - const auto end = std::chrono::system_clock::now(); + const auto end = std::chrono::system_clock::now(); - const std::chrono::duration duration = end - start; - const time_t end_time = std::chrono::system_clock::to_time_t(end); - std::cout << "End: " << ctime(&end_time); - std::cout << "Duration: " << duration.count() << "s\n" << std::endl; + const std::chrono::duration duration = end - start; + const time_t end_time = std::chrono::system_clock::to_time_t(end); + std::cout << "End: " << ctime(&end_time); + std::cout << "Duration: " << duration.count() << "s\n" << std::endl; } diff --git a/src/sims.h b/src/sims.h index 094ce07..4dda6d8 100644 --- a/src/sims.h +++ b/src/sims.h @@ -4,21 +4,21 @@ #include "dynamics/base.h" #include "experiments/base.h" -#include -#include #include #include +#include +#include -void run_simulation( - Experiment& experiment, - const std::unordered_map& parameter, - const std::unordered_map& optional, - Dynamics& dynamics, - std::mt19937_64& rng); +void run_simulation(Experiment &experiment, + const std::unordered_map ¶meter, + const std::unordered_map &optional, + Dynamics &dynamics, std::mt19937_64 &rng); -Trajectory make_trajectory(Dynamics& dynamics, double t_max, std::mt19937_64& rng); +Trajectory make_trajectory(Dynamics &dynamics, double t_max, + std::mt19937_64 &rng); -std::chrono::system_clock::time_point printStart(const std::unordered_map &optional); +std::chrono::system_clock::time_point +printStart(const std::unordered_map &optional); void printEnd(std::chrono::system_clock::time_point start); -#endif //RWSIM_SIMS_H +#endif // RWSIM_SIMS_H diff --git a/src/times/CMakeLists.txt b/src/times/CMakeLists.txt index 5aeca76..8f7e468 100644 --- a/src/times/CMakeLists.txt +++ b/src/times/CMakeLists.txt @@ -1,9 +1,5 @@ -add_library( - RWTime OBJECT - base.cpp - base.h - delta.cpp - delta.h - lognormal.cpp - lognormal.h +add_library(RWTime OBJECT + base.cpp base.h + delta.cpp delta.h + lognormal.cpp lognormal.h ) diff --git a/src/times/base.cpp b/src/times/base.cpp index 7a0408c..01e7b8a 100644 --- a/src/times/base.cpp +++ b/src/times/base.cpp @@ -2,26 +2,30 @@ #include "../utils/registry.h" namespace times { - BaseDistribution::BaseDistribution(std::string name, const double tau) : m_name(std::move(name)), m_tau(tau), m_tau_jump(tau) {} +BaseDistribution::BaseDistribution(std::string name, const double tau) + : m_name(std::move(name)), m_tau(tau), m_tau_jump(tau) {} - BaseDistribution::BaseDistribution(std::string name) : m_name(std::move(name)) {} +BaseDistribution::BaseDistribution(std::string name) + : m_name(std::move(name)) {} - double BaseDistribution::tau_wait(std::mt19937_64& rng) const { - return std::exponential_distribution(1./m_tau_jump)(rng); - } - - void BaseDistribution::setParameters(const std::unordered_map ¶meters) { - m_tau = parameters.at("tau"); - } - - std::unordered_map BaseDistribution::getParameters() const { - return std::unordered_map{ - {"tau", m_tau}, - }; - } - - - std::unique_ptr BaseDistribution::createFromInput(const std::string& input) { - return Registry::create(input); - } +double BaseDistribution::tau_wait(std::mt19937_64 &rng) const { + return std::exponential_distribution(1. / m_tau_jump)(rng); } + +void BaseDistribution::setParameters( + const std::unordered_map ¶meters) { + m_tau = parameters.at("tau"); +} + +std::unordered_map +BaseDistribution::getParameters() const { + return std::unordered_map{ + {"tau", m_tau}, + }; +} + +std::unique_ptr +BaseDistribution::createFromInput(const std::string &input) { + return Registry::create(input); +} +} // namespace times diff --git a/src/times/base.h b/src/times/base.h index 1f23c68..bbe7aa0 100644 --- a/src/times/base.h +++ b/src/times/base.h @@ -6,33 +6,35 @@ #include namespace times { - class BaseDistribution { - public: - virtual ~BaseDistribution() = default; +class BaseDistribution { +public: + virtual ~BaseDistribution() = default; - BaseDistribution(std::string, double); - explicit BaseDistribution(std::string); + BaseDistribution(std::string, double); + explicit BaseDistribution(std::string); - [[nodiscard]] double getTau() const { return m_tau; } - void setTau(const double tau) { m_tau = tau; } - [[nodiscard]] std::string getName() const { return m_name; } + [[nodiscard]] double getTau() const { return m_tau; } + void setTau(const double tau) { m_tau = tau; } + [[nodiscard]] std::string getName() const { return m_name; } - virtual void setParameters(const std::unordered_map&); - [[nodiscard]] virtual std::unordered_map getParameters() const; + virtual void setParameters(const std::unordered_map &); + [[nodiscard]] virtual std::unordered_map + getParameters() const; - virtual void initialize(std::mt19937_64& rng) = 0; - [[nodiscard]] double tau_wait(std::mt19937_64& rng) const; - [[nodiscard]] virtual std::unique_ptr clone() const = 0; + virtual void initialize(std::mt19937_64 &rng) = 0; + [[nodiscard]] double tau_wait(std::mt19937_64 &rng) const; + [[nodiscard]] virtual std::unique_ptr clone() const = 0; - [[nodiscard]] virtual std::string toString() const = 0; + [[nodiscard]] virtual std::string toString() const = 0; - static std::unique_ptr createFromInput(const std::string& input); + static std::unique_ptr + createFromInput(const std::string &input); - protected: - std::string m_name{"BaseDistribution"}; - double m_tau{1.}; - double m_tau_jump{1.}; - }; -} +protected: + std::string m_name{"BaseDistribution"}; + double m_tau{1.}; + double m_tau_jump{1.}; +}; +} // namespace times -#endif //RWSIM_TIMESBASE_H +#endif // RWSIM_TIMESBASE_H diff --git a/src/times/delta.cpp b/src/times/delta.cpp index 93efbd4..8383c77 100644 --- a/src/times/delta.cpp +++ b/src/times/delta.cpp @@ -1,21 +1,22 @@ #include "delta.h" #include "../utils/registry.h" -static AutoRegister reg("Delta"); +static AutoRegister + reg("Delta"); namespace times { - DeltaDistribution::DeltaDistribution(const double tau) : BaseDistribution(std::string("Delta"), tau) {} - DeltaDistribution::DeltaDistribution() : BaseDistribution(std::string("Delta")) {} +DeltaDistribution::DeltaDistribution(const double tau) + : BaseDistribution(std::string("Delta"), tau) {} +DeltaDistribution::DeltaDistribution() + : BaseDistribution(std::string("Delta")) {} - void DeltaDistribution::initialize(std::mt19937_64&) { - m_tau_jump = m_tau; - } +void DeltaDistribution::initialize(std::mt19937_64 &) { m_tau_jump = m_tau; } - std::unique_ptr DeltaDistribution::clone() const { - return std::make_unique(*this); - } - - std::string DeltaDistribution::toString() const { - return {"Delta/tau=" + std::to_string(m_tau)}; - } +std::unique_ptr DeltaDistribution::clone() const { + return std::make_unique(*this); } + +std::string DeltaDistribution::toString() const { + return {"Delta/tau=" + std::to_string(m_tau)}; +} +} // namespace times diff --git a/src/times/delta.h b/src/times/delta.h index 30a57b5..82fce55 100644 --- a/src/times/delta.h +++ b/src/times/delta.h @@ -4,15 +4,15 @@ #include "base.h" namespace times { - class DeltaDistribution final : public BaseDistribution { - public: - explicit DeltaDistribution(double); - DeltaDistribution(); +class DeltaDistribution final : public BaseDistribution { +public: + explicit DeltaDistribution(double); + DeltaDistribution(); - void initialize(std::mt19937_64& rng) override; - [[nodiscard]] std::unique_ptr clone() const override; + void initialize(std::mt19937_64 &rng) override; + [[nodiscard]] std::unique_ptr clone() const override; - [[nodiscard]] std::string toString() const override; - }; -} -#endif //RWSIM_TIMESDELTA_H + [[nodiscard]] std::string toString() const override; +}; +} // namespace times +#endif // RWSIM_TIMESDELTA_H diff --git a/src/times/lognormal.cpp b/src/times/lognormal.cpp index 7d64f16..7d41442 100644 --- a/src/times/lognormal.cpp +++ b/src/times/lognormal.cpp @@ -2,33 +2,41 @@ #include "../utils/registry.h" #include -static AutoRegister reg("LogNormal"); +static AutoRegister + reg("LogNormal"); namespace times { - LogNormalDistribution::LogNormalDistribution(const double tau, const double sigma) : BaseDistribution(std::string("LogNormal"), tau), m_sigma(sigma), m_distribution(std::log(tau), sigma) {} - LogNormalDistribution::LogNormalDistribution() : BaseDistribution(std::string("LogNormal")) {} +LogNormalDistribution::LogNormalDistribution(const double tau, + const double sigma) + : BaseDistribution(std::string("LogNormal"), tau), m_sigma(sigma), + m_distribution(std::log(tau), sigma) {} +LogNormalDistribution::LogNormalDistribution() + : BaseDistribution(std::string("LogNormal")) {} - void LogNormalDistribution::setParameters(const std::unordered_map ¶meters) { - m_sigma = parameters.at("sigma"); - BaseDistribution::setParameters(parameters); - } - - std::unordered_map LogNormalDistribution::getParameters() const { - auto parameter = BaseDistribution::getParameters(); - parameter["sigma"] = m_sigma; - return parameter; - } - - void LogNormalDistribution::initialize(std::mt19937_64& rng) { - m_distribution = std::lognormal_distribution(std::log(m_tau), m_sigma); - m_tau_jump = m_distribution(rng); - } - - std::unique_ptr LogNormalDistribution::clone() const { - return std::make_unique(*this); - } - - std::string LogNormalDistribution::toString() const { - return {"LogNormal/tau=" + std::to_string(m_tau) + "/sigma=" + std::to_string(m_sigma)}; - } +void LogNormalDistribution::setParameters( + const std::unordered_map ¶meters) { + m_sigma = parameters.at("sigma"); + BaseDistribution::setParameters(parameters); } + +std::unordered_map +LogNormalDistribution::getParameters() const { + auto parameter = BaseDistribution::getParameters(); + parameter["sigma"] = m_sigma; + return parameter; +} + +void LogNormalDistribution::initialize(std::mt19937_64 &rng) { + m_distribution = std::lognormal_distribution(std::log(m_tau), m_sigma); + m_tau_jump = m_distribution(rng); +} + +std::unique_ptr LogNormalDistribution::clone() const { + return std::make_unique(*this); +} + +std::string LogNormalDistribution::toString() const { + return {"LogNormal/tau=" + std::to_string(m_tau) + + "/sigma=" + std::to_string(m_sigma)}; +} +} // namespace times diff --git a/src/times/lognormal.h b/src/times/lognormal.h index 09f789b..24cbf5d 100644 --- a/src/times/lognormal.h +++ b/src/times/lognormal.h @@ -5,23 +5,24 @@ #include namespace times { - class LogNormalDistribution final : public BaseDistribution { - public: - LogNormalDistribution(double, double); - LogNormalDistribution(); +class LogNormalDistribution final : public BaseDistribution { +public: + LogNormalDistribution(double, double); + LogNormalDistribution(); - void setParameters(const std::unordered_map &) override; - [[nodiscard]] std::unordered_map getParameters() const override; + void setParameters(const std::unordered_map &) override; + [[nodiscard]] std::unordered_map + getParameters() const override; - [[nodiscard]] std::string toString() const override; - [[nodiscard]] std::unique_ptr clone() const override; + [[nodiscard]] std::string toString() const override; + [[nodiscard]] std::unique_ptr clone() const override; - void initialize(std::mt19937_64& rng) override; + void initialize(std::mt19937_64 &rng) override; - private: - double m_sigma{1}; - std::lognormal_distribution<> m_distribution; - }; -} +private: + double m_sigma{1}; + std::lognormal_distribution<> m_distribution; +}; +} // namespace times -#endif //LOGNORMAL_H +#endif // LOGNORMAL_H diff --git a/src/utils/functions.cpp b/src/utils/functions.cpp index 2ff18e5..bfa1799 100644 --- a/src/utils/functions.cpp +++ b/src/utils/functions.cpp @@ -1,51 +1,54 @@ -#include #include #include +#include - -int nearest_index(const std::vector &x_ref, const double x, int start=0) { - const int last = static_cast(x_ref.size()) - 2; - while (start < last && x > x_ref[start+1]) { - start++; - } - return start; +int nearest_index(const std::vector &x_ref, const double x, + int start = 0) { + const int last = static_cast(x_ref.size()) - 2; + while (start < last && x > x_ref[start + 1]) { + start++; + } + return start; } -double lerp(const std::vector& x_ref, const std::vector& y_ref, const double x, const int i) { - /* - * Linear interpolation between two - */ - const double x_left = x_ref[i]; - const double y_left = y_ref[i]; - const double x_right = x_ref[i+1]; - const double y_right = y_ref[i+1]; +double lerp(const std::vector &x_ref, const std::vector &y_ref, + const double x, const int i) { + /* + * Linear interpolation between two + */ + const double x_left = x_ref[i]; + const double y_left = y_ref[i]; + const double x_right = x_ref[i + 1]; + const double y_right = y_ref[i + 1]; - const double dydx = (y_right - y_left) / ( x_right - x_left ); + const double dydx = (y_right - y_left) / (x_right - x_left); - return y_left + dydx * (x - x_left); + return y_left + dydx * (x - x_left); } - std::chrono::time_point printSteps( /* - * Prints roughly every 10 seconds how many runs were done and gives a time estimation + * Prints roughly every 10 seconds how many runs were done and gives a time + * estimation */ - const std::chrono::time_point last_print_out, + const std::chrono::time_point + last_print_out, const std::chrono::time_point start, - const int total, - const int steps - ) { - const auto now = std::chrono::high_resolution_clock::now(); + const int total, const int steps) { + const auto now = std::chrono::high_resolution_clock::now(); - if (const std::chrono::duration duration = now - last_print_out; duration.count() < 10.) { - return last_print_out; - } + if (const std::chrono::duration duration = now - last_print_out; + duration.count() < 10.) { + return last_print_out; + } - const std::chrono::duration duration = now - start; - const auto passed = duration.count(); + const std::chrono::duration duration = now - start; + const auto passed = duration.count(); + std::cout << steps << " of " << total << " steps: " << passed << "s passed; ~" + << passed * static_cast(total - steps) / + static_cast(steps) + << "s remaining\n"; - std::cout << steps << " of " << total << " steps: " << passed << "s passed; ~" << passed * static_cast(total-steps) / static_cast(steps) << "s remaining\n"; - - return now; -} \ No newline at end of file + return now; +} diff --git a/src/utils/functions.h b/src/utils/functions.h index 8bd3e8d..8c6bbf5 100644 --- a/src/utils/functions.h +++ b/src/utils/functions.h @@ -1,14 +1,16 @@ #ifndef RWSIM_FUNCTIONS_H #define RWSIM_FUNCTIONS_H -#include #include +#include +int nearest_index(const std::vector &, double, int); -int nearest_index(const std::vector&, double, int); +double lerp(const std::vector &, const std::vector &, double, + int); -double lerp(const std::vector&, const std::vector&, double, int); - -std::chrono::time_point printSteps(std::chrono::time_point, std::chrono::time_point, int, int); +std::chrono::time_point +printSteps(std::chrono::time_point, + std::chrono::time_point, int, int); #endif diff --git a/src/utils/io.cpp b/src/utils/io.cpp index 3a67a1f..71b38e0 100644 --- a/src/utils/io.cpp +++ b/src/utils/io.cpp @@ -1,28 +1,28 @@ #include "io.h" -#include -#include -#include #include -#include -#include -#include -#include -#include #include +#include +#include +#include +#include +#include +#include +#include +std::pair +get_optional_parameter(std::vector::const_iterator &it) { + std::string stripped_arg; + if (it->size() > 2 && it->at(0) == '-' && it->at(1) == '-') { + stripped_arg = it->substr(2, it->size()); + } else if (it->size() > 1 && it->at(0) == '-') { + stripped_arg = it->substr(1, it->size()); + } + std::transform(stripped_arg.begin(), stripped_arg.end(), stripped_arg.begin(), + [](unsigned char c) { return std::tolower(c); }); + const auto stripped_value = std::stod(*++it, nullptr); -std::pair get_optional_parameter(std::vector::const_iterator &it) { - std::string stripped_arg; - if (it->size() > 2 && it->at(0) == '-' && it->at(1) == '-') { - stripped_arg = it->substr(2, it->size()); - } else if (it->size() > 1 && it->at(0) == '-') { - stripped_arg = it->substr(1, it->size()); - } - std::transform(stripped_arg.begin(), stripped_arg.end(), stripped_arg.begin(), [](unsigned char c) { return std::tolower(c); }); - const auto stripped_value = std::stod(*++it, nullptr); - - return std::make_pair(stripped_arg, stripped_value); + return std::make_pair(stripped_arg, stripped_value); } /** @@ -31,201 +31,197 @@ std::pair get_optional_parameter(std::vector:: * @param argv List of command-line arguments * @return Arguments */ -Arguments parse_args(const int argc, char* argv[]) { - if (argc < 3) { - throw std::runtime_error("Not enough arguments: missing parameter file"); +Arguments parse_args(const int argc, char *argv[]) { + if (argc < 3) { + throw std::runtime_error("Not enough arguments: missing parameter file"); + } + + Arguments args; + + // convert to vector to use iterator for loop + const std::vector input_args(argv + 1, argv + argc); + + for (auto it = input_args.begin(); it != input_args.end(); ++it) { + // check for optional parameter + if (it->at(0) == '-') { + + // only --spectrum and --ste are parameter that are predefined + if (*it == "--spectrum") { + args.spectrum = true; + continue; + } + + if (*it == "--ste") { + args.ste = true; + continue; + } + + // all other arguments are optional parameter + auto [option_name, option_value] = get_optional_parameter(it); + args.optional[option_name] = option_value; + continue; } - Arguments args; - - // convert to vector to use iterator for loop - const std::vector input_args(argv + 1, argv + argc); - - for (auto it = input_args.begin(); it != input_args.end(); ++it) { - // check for optional parameter - if (it->at(0) == '-') { - - // only --spectrum and --ste are parameter that are predefined - if (*it == "--spectrum") { - args.spectrum = true; - continue; - } - - if (*it == "--ste") { - args.ste = true; - continue; - } - - // all other arguments are optional parameter - auto [option_name, option_value] = get_optional_parameter(it); - args.optional[option_name] = option_value; - 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"); + // 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.parameter_file.empty() || args.motion_type.empty() || args.distribution_type.empty()) { - throw std::invalid_argument("Missing argument"); + if (args.motion_type.empty()) { + args.motion_type = *it; + continue; } - return args; + if (args.distribution_type.empty()) { + args.distribution_type = *it; + continue; + } + + throw std::invalid_argument("too many positional arguments"); + } + + if (args.parameter_file.empty() || args.motion_type.empty() || + args.distribution_type.empty()) { + throw std::invalid_argument("Missing argument"); + } + + return args; } +std::unordered_map +read_parameter(const std::filesystem::path &infile) { + if (!exists(infile)) { + std::cerr << "File " << infile << " does not exist" << std::endl; + exit(1); + } -std::unordered_map read_parameter(const std::filesystem::path& infile) { - if (!exists(infile)) { - std::cerr << "File " << infile << " does not exist" << std::endl; - exit(1); - } + std::ifstream instream(infile); - std::ifstream instream(infile); + std::unordered_map parameter; - std::unordered_map parameter; + std::string line; + std::string delim = "="; + std::string key; + std::string value; + size_t delim_pos; - std::string line; - std::string delim = "="; - std::string key; - std::string value; - size_t delim_pos; + while (std::getline(instream, line)) { + // skip comment lines starting with #, and empty lines + if (line[0] == '#' || line.length() == 1) + continue; - while (std::getline(instream, line)) { - // skip comment lines starting with #, and empty lines - if (line[0] == '#' || line.length() == 1) continue; + // strip spaces from line to have always key=value + line.erase(std::remove(line.begin(), line.end(), ' '), line.end()); - // strip spaces from line to have always key=value - line.erase(std::remove(line.begin(), line.end(), ' '), line.end()); + // split at '=' character and add to map + delim_pos = line.find('='); + key = line.substr(0, delim_pos); + value = line.substr(delim_pos + 1); + parameter[key] = std::stod(value); + } - // split at '=' character and add to map - delim_pos = line.find('='); - key = line.substr(0, delim_pos); - value = line.substr(delim_pos+1); - parameter[key] = std::stod(value); - } - - return parameter; + return parameter; } +std::string make_directory(const std::string &path_name) { + if (!std::filesystem::create_directories(path_name)) { + std::cout << "Created directory " << path_name << std::endl; + } -std::string make_directory(const std::string& path_name) { - if (!std::filesystem::create_directories(path_name)) { - std::cout << "Created directory " << path_name << std::endl; - } - - return path_name; + return path_name; } - void save_parameter_to_file( - const std::string& resulttype, - const std::string& directory, - const std::unordered_map& parameter, - const std::unordered_map& optional - ) { - std::ostringstream parameter_filename; - parameter_filename << directory << "/" << resulttype; + const std::string &resulttype, const std::string &directory, + const std::unordered_map ¶meter, + const std::unordered_map &optional) { + std::ostringstream parameter_filename; + parameter_filename << directory << "/" << resulttype; - parameter_filename << std::setprecision(6) << std::scientific; - for (const auto& [key, value]: optional) { - parameter_filename << "_" << key << "=" << value; - } - parameter_filename << "_parameter.txt"; - - { - // write data to file, columns are secondary axis (echo delay, evolution times) - std::string datafile = parameter_filename.str(); - std::ofstream filestream(datafile, std::ios::out); - - for (const auto& [key, value]: parameter) { - filestream << key << "=" << value << "\n"; - } - } -} - - -void save_data_to_file( - const std::string& resulttype, - const std::string& directory, - const std::vector& x, - const std::map>& y, - const std::unordered_map& optional - ) { - // make file name - std::ostringstream datafile_name; - datafile_name << directory << "/" << resulttype; - - datafile_name << std::setprecision(6) << std::scientific; - for (const auto& [key, value]: optional) { - datafile_name << "_" << key << "=" << value; - } - datafile_name << ".dat"; - - { - // write data to file, columns are secondary axis (echo delay, evolution times) - std::string datafile = datafile_name.str(); - std::ofstream filestream(datafile, std::ios::out); - - // first line are values of secondary axis - filestream << "#"; - for (const auto& [t_echo_j, _] : y) { - filestream << t_echo_j << "\t"; - } - filestream << std::endl; - - // write values to file - auto size = x.size(); - for (unsigned int i = 0; i < size; i++) { - filestream << x[i]; - for (const auto& [_, fid_j] : y) { - filestream << "\t" << fid_j[i]; - } - filestream << "\n"; - } + parameter_filename << std::setprecision(6) << std::scientific; + for (const auto &[key, value] : optional) { + parameter_filename << "_" << key << "=" << value; + } + parameter_filename << "_parameter.txt"; + + { + // write data to file, columns are secondary axis (echo delay, evolution + // times) + std::string datafile = parameter_filename.str(); + std::ofstream filestream(datafile, std::ios::out); + + for (const auto &[key, value] : parameter) { + filestream << key << "=" << value << "\n"; } + } } void save_data_to_file( - const std::string& resulttype, - const std::string& directory, - const std::vector& x, - const std::vector& y, - const std::unordered_map& optional - ) { - // make file name - std::ostringstream datafile_name; - datafile_name << directory << "/" << resulttype; + const std::string &resulttype, const std::string &directory, + const std::vector &x, + const std::map> &y, + const std::unordered_map &optional) { + // make file name + std::ostringstream datafile_name; + datafile_name << directory << "/" << resulttype; - datafile_name << std::setprecision(6) << std::scientific; - for (const auto& [key, value]: optional) { - datafile_name << "_" << key << "=" << value; + datafile_name << std::setprecision(6) << std::scientific; + for (const auto &[key, value] : optional) { + datafile_name << "_" << key << "=" << value; + } + datafile_name << ".dat"; + + { + // write data to file, columns are secondary axis (echo delay, evolution + // times) + std::string datafile = datafile_name.str(); + std::ofstream filestream(datafile, std::ios::out); + + // first line are values of secondary axis + filestream << "#"; + for (const auto &[t_echo_j, _] : y) { + filestream << t_echo_j << "\t"; } - datafile_name << ".dat"; + filestream << std::endl; - { - // write data to file, columns are secondary axis (echo delay, evolution times) - std::string datafile = datafile_name.str(); - std::ofstream filestream(datafile, std::ios::out); - - // write values to file - auto size = x.size(); - for (unsigned int i = 0; i < size; i++) { - filestream << x[i] << "\t" << y[i] << "\n"; - } + // write values to file + auto size = x.size(); + for (unsigned int i = 0; i < size; i++) { + filestream << x[i]; + for (const auto &[_, fid_j] : y) { + filestream << "\t" << fid_j[i]; + } + filestream << "\n"; } + } +} + +void save_data_to_file( + const std::string &resulttype, const std::string &directory, + const std::vector &x, const std::vector &y, + const std::unordered_map &optional) { + // make file name + std::ostringstream datafile_name; + datafile_name << directory << "/" << resulttype; + + datafile_name << std::setprecision(6) << std::scientific; + for (const auto &[key, value] : optional) { + datafile_name << "_" << key << "=" << value; + } + datafile_name << ".dat"; + + { + // write data to file, columns are secondary axis (echo delay, evolution + // times) + std::string datafile = datafile_name.str(); + std::ofstream filestream(datafile, std::ios::out); + + // write values to file + auto size = x.size(); + for (unsigned int i = 0; i < size; i++) { + filestream << x[i] << "\t" << y[i] << "\n"; + } + } } diff --git a/src/utils/io.h b/src/utils/io.h index 50ab2e4..e0cadf0 100644 --- a/src/utils/io.h +++ b/src/utils/io.h @@ -1,29 +1,38 @@ #ifndef RWSIM_IO_H #define RWSIM_IO_H -#include +#include #include #include -#include +#include #include struct Arguments { - std::string parameter_file{}; - bool ste = false; - bool spectrum = false; - std::string motion_type{}; - std::string distribution_type{}; - std::unordered_map optional; + std::string parameter_file{}; + bool ste = false; + bool spectrum = false; + std::string motion_type{}; + std::string distribution_type{}; + std::unordered_map optional; }; -Arguments parse_args(int argc, char* argv[]); -std::pair get_optional_parameter(std::vector::const_iterator &it); -std::unordered_map read_parameter(const std::filesystem::path&); +Arguments parse_args(int argc, char *argv[]); +std::pair +get_optional_parameter(std::vector::const_iterator &it); +std::unordered_map +read_parameter(const std::filesystem::path &); -std::string make_directory(const std::string& path_name); +std::string make_directory(const std::string &path_name); -void save_parameter_to_file(const std::string&, const std::string&, const std::unordered_map&, const std::unordered_map&); -void save_data_to_file(const std::string&, const std::string&, const std::vector&, const std::map>&, const std::unordered_map&); -void save_data_to_file(const std::string&, const std::string&, const std::vector&, const std::vector&, const std::unordered_map&); +void save_parameter_to_file(const std::string &, const std::string &, + const std::unordered_map &, + const std::unordered_map &); +void save_data_to_file(const std::string &, const std::string &, + const std::vector &, + const std::map> &, + const std::unordered_map &); +void save_data_to_file(const std::string &, const std::string &, + const std::vector &, const std::vector &, + const std::unordered_map &); #endif diff --git a/src/utils/ranges.cpp b/src/utils/ranges.cpp index 732c4f0..cf0ad95 100644 --- a/src/utils/ranges.cpp +++ b/src/utils/ranges.cpp @@ -1,57 +1,58 @@ -#include #include #include +#include #include "ranges.h" +std::vector arange(const int size, const double spacing = 1.) { + std::vector out(size); + std::generate(out.begin(), out.end(), + [n = 0, spacing]() mutable { return n++ * spacing; }); -std::vector arange(const int size, const double spacing=1.) { - std::vector out(size); - std::generate(out.begin(), out.end(), [n = 0, spacing]() mutable { return n++ * spacing; }); - - return out; + return out; } -std::vector linspace(const double start, const double stop, const int steps) { - std::vector range; - range.reserve(steps); +std::vector linspace(const double start, const double stop, + const int steps) { + std::vector range; + range.reserve(steps); - if (steps == 0) { - return range; - } - - if (steps == 1) { - range.push_back(start); - return range; - } - - const double stepsize = (stop-start) / (steps-1); - for (int i=0; i logspace(const double start, const double stop, + const int steps) { + std::vector range; + range.reserve(steps); -std::vector logspace(const double start, const double stop, const int steps) { - std::vector range; - range.reserve(steps); - - if (steps == 0) { - return range; - } - - if (steps == 1) { - range.push_back(start); - return range; - } - - const double logstart = std::log10(start); - const double logstop = std::log10(stop); - - const double stepsize = (logstop-logstart) / (steps-1); - for (int i=0; i #include +#include #include #include -#include #include -template -class Registry { +template class Registry { public: - using Creator = std::function()>; + using Creator = std::function()>; - static std::unordered_map& entries() { - static std::unordered_map map; - return map; - } + static std::unordered_map &entries() { + static std::unordered_map map; + return map; + } - static void add(const std::string& name, Creator creator) { - entries()[name] = std::move(creator); - } + static void add(const std::string &name, Creator creator) { + entries()[name] = std::move(creator); + } - static std::unique_ptr create(const std::string& name) { - auto& map = entries(); - auto it = map.find(name); - if (it == map.end()) { - std::string msg = "Unknown model '" + name + "'. Available: "; - for (const auto& [key, _] : map) { - msg += key + ", "; - } - if (!map.empty()) { - msg.pop_back(); - msg.pop_back(); - } - throw std::invalid_argument(msg); - } - return it->second(); + static std::unique_ptr create(const std::string &name) { + auto &map = entries(); + auto it = map.find(name); + if (it == map.end()) { + std::string msg = "Unknown model '" + name + "'. Available: "; + for (const auto &[key, _] : map) { + msg += key + ", "; + } + if (!map.empty()) { + msg.pop_back(); + msg.pop_back(); + } + throw std::invalid_argument(msg); } + return it->second(); + } }; -template -struct AutoRegister { - explicit AutoRegister(const std::string& name) { - Registry::add(name, []() { return std::make_unique(); }); - } +template struct AutoRegister { + explicit AutoRegister(const std::string &name) { + Registry::add(name, []() { return std::make_unique(); }); + } }; -#endif //RWSIM_REGISTRY_H +#endif // RWSIM_REGISTRY_H