NMRRWSims/cpp
1
0
Fork 0
Code Issues 3 Pull Requests Actions Packages Projects Releases Wiki Activity

Formatting

Browse Source
This commit is contained in:
Dominik Demuth
2026-03-08 14:01:37 +01:00
parent 6579bb028b
commit 285c78bed5
52 changed files with 1381 additions and 1328 deletions
Download Patch File Download Diff File Expand all files Collapse all files

26
src/CMakeLists.txt
View File

@@ -5,14 +5,24 @@ add_subdirectory(utils)
add_subdirectory(experiments) add_subdirectory(experiments)
add_subdirectory(dynamics) add_subdirectory(dynamics)
add_library(simulation STATIC sims.cpp sims.h) add_library(simulation STATIC
target_link_libraries(simulation PRIVATE utils experiments dynamics OpenMP::OpenMP_CXX) sims.cpp sims.h
)
add_executable( target_link_libraries(simulation PRIVATE
rwsim utils
main.cpp 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)

8
src/dynamics/CMakeLists.txt
View File

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

21
src/dynamics/base.h
View File

@@ -7,20 +7,21 @@
#include <unordered_map> #include <unordered_map>
struct Step { struct Step {
double dt; double dt;
double omega; double omega;
}; };
class Dynamics { class Dynamics {
public: public:
virtual ~Dynamics() = default; virtual ~Dynamics() = default;
virtual void setParameters(const std::unordered_map<std::string, double>& parameters) = 0; virtual void
virtual void initialize(std::mt19937_64& rng) = 0; setParameters(const std::unordered_map<std::string, double> &parameters) = 0;
virtual Step next(std::mt19937_64& rng) = 0; virtual void initialize(std::mt19937_64 &rng) = 0;
[[nodiscard]] virtual double getInitOmega() const = 0; virtual Step next(std::mt19937_64 &rng) = 0;
[[nodiscard]] virtual std::unique_ptr<Dynamics> clone() const = 0; [[nodiscard]] virtual double getInitOmega() const = 0;
[[nodiscard]] virtual std::string toString() const = 0; [[nodiscard]] virtual std::unique_ptr<Dynamics> clone() const = 0;
[[nodiscard]] virtual std::string toString() const = 0;
}; };
#endif //RWSIM_DYNAMICSBASE_H #endif // RWSIM_DYNAMICSBASE_H

29
src/dynamics/decoupled.cpp
View File

@@ -5,31 +5,32 @@ DecoupledDynamics::DecoupledDynamics(
std::unique_ptr<times::BaseDistribution> dist) std::unique_ptr<times::BaseDistribution> dist)
: m_motion(std::move(motion)), m_dist(std::move(dist)) {} : m_motion(std::move(motion)), m_dist(std::move(dist)) {}
void DecoupledDynamics::setParameters(const std::unordered_map<std::string, double>& parameters) { void DecoupledDynamics::setParameters(
m_motion->setParameters(parameters); const std::unordered_map<std::string, double> &parameters) {
m_dist->setParameters(parameters); m_motion->setParameters(parameters);
m_dist->setParameters(parameters);
} }
void DecoupledDynamics::initialize(std::mt19937_64& rng) { void DecoupledDynamics::initialize(std::mt19937_64 &rng) {
m_motion->initialize(rng); m_motion->initialize(rng);
m_dist->initialize(rng); m_dist->initialize(rng);
} }
Step DecoupledDynamics::next(std::mt19937_64& rng) { Step DecoupledDynamics::next(std::mt19937_64 &rng) {
double dt = m_dist->tau_wait(rng); double dt = m_dist->tau_wait(rng);
double omega = m_motion->jump(rng); double omega = m_motion->jump(rng);
return {dt, omega}; return {dt, omega};
} }
double DecoupledDynamics::getInitOmega() const { double DecoupledDynamics::getInitOmega() const {
return m_motion->getInitOmega(); return m_motion->getInitOmega();
} }
std::unique_ptr<Dynamics> DecoupledDynamics::clone() const { std::unique_ptr<Dynamics> DecoupledDynamics::clone() const {
return std::make_unique<DecoupledDynamics>( return std::make_unique<DecoupledDynamics>(m_motion->clone(),
m_motion->clone(), m_dist->clone()); m_dist->clone());
} }
std::string DecoupledDynamics::toString() const { std::string DecoupledDynamics::toString() const {
return m_motion->toString() + "/" + m_dist->toString(); return m_motion->toString() + "/" + m_dist->toString();
} }

25
src/dynamics/decoupled.h
View File

@@ -1,27 +1,28 @@
#ifndef RWSIM_DECOUPLED_H #ifndef RWSIM_DECOUPLED_H
#define RWSIM_DECOUPLED_H #define RWSIM_DECOUPLED_H
#include "base.h"
#include "../motions/base.h" #include "../motions/base.h"
#include "../times/base.h" #include "../times/base.h"
#include "base.h"
#include <memory> #include <memory>
class DecoupledDynamics final : public Dynamics { class DecoupledDynamics final : public Dynamics {
public: public:
DecoupledDynamics(std::unique_ptr<motions::BaseMotion> motion, DecoupledDynamics(std::unique_ptr<motions::BaseMotion> motion,
std::unique_ptr<times::BaseDistribution> dist); std::unique_ptr<times::BaseDistribution> dist);
void setParameters(const std::unordered_map<std::string, double>& parameters) override; void setParameters(
void initialize(std::mt19937_64& rng) override; const std::unordered_map<std::string, double> &parameters) override;
Step next(std::mt19937_64& rng) override; void initialize(std::mt19937_64 &rng) override;
[[nodiscard]] double getInitOmega() const override; Step next(std::mt19937_64 &rng) override;
[[nodiscard]] std::unique_ptr<Dynamics> clone() const override; [[nodiscard]] double getInitOmega() const override;
[[nodiscard]] std::string toString() const override; [[nodiscard]] std::unique_ptr<Dynamics> clone() const override;
[[nodiscard]] std::string toString() const override;
private: private:
std::unique_ptr<motions::BaseMotion> m_motion; std::unique_ptr<motions::BaseMotion> m_motion;
std::unique_ptr<times::BaseDistribution> m_dist; std::unique_ptr<times::BaseDistribution> m_dist;
}; };
#endif //RWSIM_DECOUPLED_H #endif // RWSIM_DECOUPLED_H

11
src/experiments/CMakeLists.txt
View File

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

28
src/experiments/base.h
View File

@@ -2,27 +2,29 @@
#define RWSIM_EXPERIMENTBASE_H #define RWSIM_EXPERIMENTBASE_H
#include <memory> #include <memory>
#include <unordered_map>
#include <string> #include <string>
#include <unordered_map>
#include <vector> #include <vector>
struct Trajectory { struct Trajectory {
std::vector<double> time; std::vector<double> time;
std::vector<double> phase; std::vector<double> phase;
std::vector<double> omega; std::vector<double> omega;
}; };
class Experiment { class Experiment {
public: public:
virtual ~Experiment() = default; virtual ~Experiment() = default;
virtual void setup(const std::unordered_map<std::string, double>& parameter, virtual void
const std::unordered_map<std::string, double>& optional) = 0; setup(const std::unordered_map<std::string, double> &parameter,
[[nodiscard]] virtual double tmax() const = 0; const std::unordered_map<std::string, double> &optional) = 0;
virtual void accumulate(const Trajectory& traj, double init_omega, int num_walkers) = 0; [[nodiscard]] virtual double tmax() const = 0;
virtual void save(const std::string& directory) = 0; virtual void accumulate(const Trajectory &traj, double init_omega,
[[nodiscard]] virtual std::unique_ptr<Experiment> clone() const = 0; int num_walkers) = 0;
virtual void merge(const Experiment& other) = 0; virtual void save(const std::string &directory) = 0;
[[nodiscard]] virtual std::unique_ptr<Experiment> clone() const = 0;
virtual void merge(const Experiment &other) = 0;
}; };
#endif //RWSIM_EXPERIMENTBASE_H #endif // RWSIM_EXPERIMENTBASE_H

96
src/experiments/spectrum.cpp
View File

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

33
src/experiments/spectrum.h
View File

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

164
src/experiments/ste.cpp
View File

@@ -1,98 +1,108 @@
#include "ste.h" #include "ste.h"
#include "../utils/functions.h" #include "../utils/functions.h"
#include "../utils/ranges.h"
#include "../utils/io.h" #include "../utils/io.h"
#include "../utils/ranges.h"
#include <algorithm> #include <algorithm>
#include <cmath> #include <cmath>
void STEExperiment::setup( void STEExperiment::setup(
const std::unordered_map<std::string, double>& parameter, const std::unordered_map<std::string, double> &parameter,
const std::unordered_map<std::string, double>& optional const std::unordered_map<std::string, double> &optional) {
) { m_parameter = parameter;
m_parameter = parameter; m_optional = optional;
m_optional = optional; m_num_mix_times = static_cast<int>(parameter.at("tmix_steps"));
m_num_mix_times = static_cast<int>(parameter.at("tmix_steps")); m_evolution_times =
m_evolution_times = linspace(parameter.at("tevo_start"), parameter.at("tevo_stop"), static_cast<int>(parameter.at("tevo_steps"))); linspace(parameter.at("tevo_start"), parameter.at("tevo_stop"),
m_mixing_times = logspace(parameter.at("tmix_start"), parameter.at("tmix_stop"), m_num_mix_times); static_cast<int>(parameter.at("tevo_steps")));
m_tpulse4 = parameter.at("tpulse4"); 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_cc_dict.clear();
m_ss_dict.clear(); m_ss_dict.clear();
for (auto t_evo_i: m_evolution_times) { for (auto t_evo_i : m_evolution_times) {
m_cc_dict[t_evo_i] = std::vector<double>(m_num_mix_times, 0.); m_cc_dict[t_evo_i] = std::vector<double>(m_num_mix_times, 0.);
m_ss_dict[t_evo_i] = std::vector<double>(m_num_mix_times, 0.); m_ss_dict[t_evo_i] = std::vector<double>(m_num_mix_times, 0.);
} }
m_f2.assign(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 m_tmax =
+ *std::max_element(m_mixing_times.begin(), m_mixing_times.end()) *std::max_element(m_evolution_times.begin(), m_evolution_times.end()) *
+ 2 * m_tpulse4; 2 +
*std::max_element(m_mixing_times.begin(), m_mixing_times.end()) +
2 * m_tpulse4;
} }
double STEExperiment::tmax() const { double STEExperiment::tmax() const { return m_tmax; }
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) { void STEExperiment::save(const std::string &directory) {
int f2_pos = 0; save_parameter_to_file(std::string("ste"), directory, m_parameter,
for (int f2_idx = 0; f2_idx < m_num_mix_times; f2_idx++) { m_optional);
const double t_mix_f2 = m_mixing_times[f2_idx]; save_data_to_file(std::string("coscos"), directory, m_mixing_times, m_cc_dict,
f2_pos = nearest_index(traj.time, t_mix_f2, f2_pos); m_optional);
m_f2[f2_idx] += traj.omega[f2_pos] * init_omega / num_walkers; 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,
for (auto& [t_evo_j, _] : m_cc_dict) { m_optional);
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);
} }
std::unique_ptr<Experiment> STEExperiment::clone() const { std::unique_ptr<Experiment> STEExperiment::clone() const {
return std::make_unique<STEExperiment>(*this); return std::make_unique<STEExperiment>(*this);
} }
void STEExperiment::merge(const Experiment& other) { void STEExperiment::merge(const Experiment &other) {
const auto& o = dynamic_cast<const STEExperiment&>(other); const auto &o = dynamic_cast<const STEExperiment &>(other);
for (auto& [t_evo, cc] : m_cc_dict) { for (auto &[t_evo, cc] : m_cc_dict) {
const auto& other_cc = o.m_cc_dict.at(t_evo); const auto &other_cc = o.m_cc_dict.at(t_evo);
const auto& other_ss = o.m_ss_dict.at(t_evo); const auto &other_ss = o.m_ss_dict.at(t_evo);
auto& ss = m_ss_dict[t_evo]; auto &ss = m_ss_dict[t_evo];
for (size_t i = 0; i < cc.size(); i++) { for (size_t i = 0; i < cc.size(); i++) {
cc[i] += other_cc[i]; cc[i] += other_cc[i];
ss[i] += other_ss[i]; ss[i] += other_ss[i];
}
}
for (size_t i = 0; i < m_f2.size(); i++) {
m_f2[i] += o.m_f2[i];
} }
}
for (size_t i = 0; i < m_f2.size(); i++) {
m_f2[i] += o.m_f2[i];
}
} }

39
src/experiments/ste.h
View File

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

83
src/main.cpp
View File

@@ -1,61 +1,58 @@
#include "sims.h"
#include "dynamics/decoupled.h" #include "dynamics/decoupled.h"
#include "experiments/spectrum.h" #include "experiments/spectrum.h"
#include "experiments/ste.h" #include "experiments/ste.h"
#include "utils/io.h"
#include "motions/base.h" #include "motions/base.h"
#include "sims.h"
#include "times/base.h" #include "times/base.h"
#include "utils/io.h"
#include <iostream> #include <iostream>
#include <unordered_map>
#include <random> #include <random>
#include <unordered_map>
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[]) { // print parameter of simulation to inform user
Arguments args; std::cout << "Found parameter\n";
try { for (const auto &[key, value] : parameter) {
args = parse_args(argc, argv); std::cout << key << ": " << std::to_string(value) << "\n";
} catch (std::invalid_argument& error) { }
std::cerr << error.what() << std::endl; std::cout << std::endl;
return 1;
}
std::unordered_map parameter { read_parameter(args.parameter_file) }; std::mt19937_64 rng;
if (parameter.count("seed")) {
rng.seed(static_cast<uint64_t>(parameter.at("seed")));
} else {
std::random_device rd;
rng.seed(rd());
}
for (const auto& [key, value]: args.optional) { DecoupledDynamics dynamics(
parameter[key] = value; motions::BaseMotion::createFromInput(args.motion_type),
} times::BaseDistribution::createFromInput(args.distribution_type));
// print parameter of simulation to inform user if (args.spectrum) {
std::cout << "Found parameter\n"; SpectrumExperiment experiment;
for (const auto& [key, value]: parameter) { run_simulation(experiment, parameter, args.optional, dynamics, rng);
std::cout << key << ": " << std::to_string(value) << "\n"; }
} if (args.ste) {
std::cout << std::endl; STEExperiment experiment;
run_simulation(experiment, parameter, args.optional, dynamics, rng);
}
std::mt19937_64 rng; return 0;
if (parameter.count("seed")) {
rng.seed(static_cast<uint64_t>(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;
} }

40
src/motions/CMakeLists.txt
View File

@@ -1,29 +1,15 @@
# Create a library target for motions add_library(RWMotion OBJECT
add_library( base.cpp base.h
RWMotion OBJECT bimodalangle.cpp bimodalangle.h
conewobble.cpp conewobble.cpp conewobble.h
conewobble.h conemotion.cpp conemotion.h
conemotion.cpp coordinates.cpp coordinates.h
conemotion.h diffusivemotion.cpp diffusivemotion.h
diffusivemotion.cpp foursitejump.cpp foursitejump.h
diffusivemotion.h isosmallangle.cpp isosmallangle.h
coordinates.cpp nsiteconejump.cpp nsiteconejump.h
coordinates.h random.cpp random.h
base.cpp rjoac.cpp rjoac.h
base.h sixsitejump.cpp sixsitejump.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
) )

85
src/motions/base.cpp
View File

@@ -5,48 +5,51 @@
#include <iostream> #include <iostream>
namespace motions { 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 { double BaseMotion::omega_q(const double cos_theta, const double phi) const {
const double cos_theta_square = cos_theta * cos_theta; const double cos_theta_square = cos_theta * cos_theta;
const double sin_theta_square = 1. - cos_theta_square; 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)); 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> BaseMotion::createFromInput(const std::string& input) {
return Registry<BaseMotion>::create(input);
}
void BaseMotion::setParameters(const std::unordered_map<std::string, double> &parameters) {
m_delta = parameters.at("delta");
m_eta = parameters.at("eta");
}
std::unordered_map<std::string, double> BaseMotion::getParameters() const {
return std::unordered_map<std::string, double>{
{"delta", m_delta},
{"eta", m_eta}
};
}
std::ostream& operator<<(std::ostream& os, const BaseMotion& m) {
os << m.getName();
return os;
}
} }
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>
BaseMotion::createFromInput(const std::string &input) {
return Registry<BaseMotion>::create(input);
}
void BaseMotion::setParameters(
const std::unordered_map<std::string, double> &parameters) {
m_delta = parameters.at("delta");
m_eta = parameters.at("eta");
}
std::unordered_map<std::string, double> BaseMotion::getParameters() const {
return std::unordered_map<std::string, double>{{"delta", m_delta},
{"eta", m_eta}};
}
std::ostream &operator<<(std::ostream &os, const BaseMotion &m) {
os << m.getName();
return os;
}
} // namespace motions

63
src/motions/base.h
View File

@@ -8,42 +8,43 @@
#include <unordered_map> #include <unordered_map>
namespace motions { namespace motions {
class BaseMotion { class BaseMotion {
public: public:
virtual ~BaseMotion() = default; virtual ~BaseMotion() = default;
BaseMotion(std::string, double, double); BaseMotion(std::string, double, double);
explicit BaseMotion(std::string); explicit BaseMotion(std::string);
coordinates::SphericalPos draw_position(std::mt19937_64& rng); coordinates::SphericalPos draw_position(std::mt19937_64 &rng);
[[nodiscard]] double omega_q(double, double) const; [[nodiscard]] double omega_q(double, double) const;
[[nodiscard]] double omega_q(const coordinates::SphericalPos&) const; [[nodiscard]] double omega_q(const coordinates::SphericalPos &) const;
virtual void initialize(std::mt19937_64& rng) = 0; virtual void initialize(std::mt19937_64 &rng) = 0;
virtual double jump(std::mt19937_64& rng) = 0; virtual double jump(std::mt19937_64 &rng) = 0;
[[nodiscard]] virtual std::unique_ptr<BaseMotion> clone() const = 0; [[nodiscard]] virtual std::unique_ptr<BaseMotion> clone() const = 0;
virtual void setParameters(const std::unordered_map<std::string, double>&); virtual void setParameters(const std::unordered_map<std::string, double> &);
[[nodiscard]] virtual std::unordered_map<std::string, double> getParameters() const; [[nodiscard]] virtual std::unordered_map<std::string, double>
[[nodiscard]] double getDelta() const { return m_delta; } getParameters() const;
void setDelta(const double delta) { m_delta = delta; } [[nodiscard]] double getDelta() const { return m_delta; }
[[nodiscard]] double getEta() const { return m_eta; } void setDelta(const double delta) { m_delta = delta; }
void setEta(const double eta) { m_eta = eta; } [[nodiscard]] double getEta() const { return m_eta; }
[[nodiscard]] std::string getName() const { return m_name; } void setEta(const double eta) { m_eta = eta; }
[[nodiscard]] double getInitOmega() const { return m_initial_omega; } [[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<BaseMotion> createFromInput(const std::string& input); static std::unique_ptr<BaseMotion> createFromInput(const std::string &input);
protected: protected:
std::string m_name{"BaseMotion"}; std::string m_name{"BaseMotion"};
double m_delta{1.}; double m_delta{1.};
double m_eta{0.}; double m_eta{0.};
std::uniform_real_distribution<> m_uni_dist{0., 1.}; std::uniform_real_distribution<> m_uni_dist{0., 1.};
double m_initial_omega{0.}; double m_initial_omega{0.};
}; };
std::ostream& operator<<(std::ostream& os, const BaseMotion& m); std::ostream &operator<<(std::ostream &os, const BaseMotion &m);
} } // namespace motions
#endif //RWSIM_MOTIONBASE_H #endif // RWSIM_MOTIONBASE_H

82
src/motions/bimodalangle.cpp
View File

@@ -1,47 +1,53 @@
#include "bimodalangle.h" #include "bimodalangle.h"
#include "coordinates.h"
#include "../utils/registry.h" #include "../utils/registry.h"
#include "coordinates.h"
static AutoRegister<motions::BaseMotion, motions::BimodalAngle> reg("BimodalAngle"); static AutoRegister<motions::BaseMotion, motions::BimodalAngle>
reg("BimodalAngle");
namespace motions { namespace motions {
BimodalAngle::BimodalAngle(const double delta, const double eta, const double angle1, const double angle2, const double prob) : BimodalAngle::BimodalAngle(const double delta, const double eta,
DiffusiveMotion(std::string("BimodalAngle"), delta, eta), const double angle1, const double angle2,
m_angle1(angle1 * M_PI / 180.0), const double prob)
m_angle2(angle2 * M_PI / 180.0), : DiffusiveMotion(std::string("BimodalAngle"), delta, eta),
m_prob(prob) {} m_angle1(angle1 * M_PI / 180.0), m_angle2(angle2 * M_PI / 180.0),
BimodalAngle::BimodalAngle() : DiffusiveMotion(std::string("BimodalAngle")) {} m_prob(prob) {}
BimodalAngle::BimodalAngle() : DiffusiveMotion(std::string("BimodalAngle")) {}
double BimodalAngle::jump(std::mt19937_64& rng) { double BimodalAngle::jump(std::mt19937_64 &rng) {
const double angle = m_uni_dist(rng) < m_prob ? m_angle1 : m_angle2; const double angle = m_uni_dist(rng) < m_prob ? m_angle1 : m_angle2;
const double gamma{2 * M_PI * m_uni_dist(rng)}; const double gamma{2 * M_PI * m_uni_dist(rng)};
m_prev_pos = rotate(m_prev_pos, angle, gamma); m_prev_pos = rotate(m_prev_pos, angle, gamma);
return omega_q(m_prev_pos); return omega_q(m_prev_pos);
}
std::unique_ptr<BaseMotion> BimodalAngle::clone() const {
return std::make_unique<BimodalAngle>(*this);
}
void BimodalAngle::setParameters(const std::unordered_map<std::string, double> &parameter) {
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<std::string, double> 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)};
}
} }
std::unique_ptr<BaseMotion> BimodalAngle::clone() const {
return std::make_unique<BimodalAngle>(*this);
}
void BimodalAngle::setParameters(
const std::unordered_map<std::string, double> &parameter) {
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<std::string, double> 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

33
src/motions/bimodalangle.h
View File

@@ -5,21 +5,22 @@
#include "diffusivemotion.h" #include "diffusivemotion.h"
namespace motions { namespace motions {
class BimodalAngle final : public DiffusiveMotion { class BimodalAngle final : public DiffusiveMotion {
public: public:
BimodalAngle(double, double, double, double, double); BimodalAngle(double, double, double, double, double);
BimodalAngle(); BimodalAngle();
double jump(std::mt19937_64& rng) override; double jump(std::mt19937_64 &rng) override;
[[nodiscard]] std::unique_ptr<BaseMotion> clone() const override; [[nodiscard]] std::unique_ptr<BaseMotion> clone() const override;
void setParameters(const std::unordered_map<std::string, double> &) override; void setParameters(const std::unordered_map<std::string, double> &) override;
[[nodiscard]] std::unordered_map<std::string, double> getParameters() const override; [[nodiscard]] std::unordered_map<std::string, double>
[[nodiscard]] std::string toString() const override; getParameters() const override;
[[nodiscard]] std::string toString() const override;
protected: protected:
double m_angle1{0}; double m_angle1{0};
double m_angle2{0}; double m_angle2{0};
double m_prob{0}; double m_prob{0};
}; };
} } // namespace motions
#endif //BIMODALANGLE_H #endif // BIMODALANGLE_H

29
src/motions/conemotion.cpp
View File

@@ -1,18 +1,19 @@
#include "conemotion.h" #include "conemotion.h"
namespace motions { namespace motions {
void ConeMotion::initialize(std::mt19937_64& rng) { void ConeMotion::initialize(std::mt19937_64 &rng) {
m_axis = draw_position(rng); m_axis = draw_position(rng);
}
void ConeMotion::setParameters(const std::unordered_map<std::string, double> &parameters) {
BaseMotion::setParameters(parameters);
m_angle = parameters.at("angle");
}
std::unordered_map<std::string, double> ConeMotion::getParameters() const {
auto parameter = BaseMotion::getParameters();
parameter["angle"] = m_angle;
return parameter;
}
} }
void ConeMotion::setParameters(
const std::unordered_map<std::string, double> &parameters) {
BaseMotion::setParameters(parameters);
m_angle = parameters.at("angle");
}
std::unordered_map<std::string, double> ConeMotion::getParameters() const {
auto parameter = BaseMotion::getParameters();
parameter["angle"] = m_angle;
return parameter;
}
} // namespace motions

25
src/motions/conemotion.h
View File

@@ -5,19 +5,20 @@
#include "coordinates.h" #include "coordinates.h"
namespace motions { namespace motions {
class ConeMotion : public BaseMotion { class ConeMotion : public BaseMotion {
public: public:
using BaseMotion::BaseMotion; using BaseMotion::BaseMotion;
void initialize(std::mt19937_64& rng) override; void initialize(std::mt19937_64 &rng) override;
void setParameters(const std::unordered_map<std::string, double> &) override; void setParameters(const std::unordered_map<std::string, double> &) override;
[[nodiscard]] std::unordered_map<std::string, double> getParameters() const override; [[nodiscard]] std::unordered_map<std::string, double>
getParameters() const override;
protected: protected:
double m_angle{0}; double m_angle{0};
coordinates::SphericalPos m_axis{1, 0}; coordinates::SphericalPos m_axis{1, 0};
}; };
} } // namespace motions
#endif //CONEMOTION_H #endif // CONEMOTION_H

37
src/motions/conewobble.cpp
View File

@@ -1,25 +1,28 @@
#include "conewobble.h" #include "conewobble.h"
#include "coordinates.h"
#include "../utils/registry.h" #include "../utils/registry.h"
#include "coordinates.h"
static AutoRegister<motions::BaseMotion, motions::WobbleCone> reg("ConeWobble"); static AutoRegister<motions::BaseMotion, motions::WobbleCone> reg("ConeWobble");
namespace motions { namespace motions {
WobbleCone::WobbleCone(const double delta, const double eta, const double chi) : ConeMotion("Wobble in Cone", delta, eta) { m_angle = chi; } WobbleCone::WobbleCone(const double delta, const double eta, const double chi)
WobbleCone::WobbleCone() : ConeMotion("Wobble in Cone") {} : ConeMotion("Wobble in Cone", delta, eta) {
m_angle = chi;
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<BaseMotion> WobbleCone::clone() const {
return std::make_unique<WobbleCone>(*this);
}
std::string WobbleCone::toString() const {
return std::string("ConeWobble/angle=") + std::to_string(m_angle);
}
} }
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<BaseMotion> WobbleCone::clone() const {
return std::make_unique<WobbleCone>(*this);
}
std::string WobbleCone::toString() const {
return std::string("ConeWobble/angle=") + std::to_string(m_angle);
}
} // namespace motions

20
src/motions/conewobble.h
View File

@@ -4,14 +4,14 @@
#include "conemotion.h" #include "conemotion.h"
namespace motions { namespace motions {
class WobbleCone final: public ConeMotion { class WobbleCone final : public ConeMotion {
public: public:
WobbleCone(double, double, double); WobbleCone(double, double, double);
WobbleCone(); WobbleCone();
double jump(std::mt19937_64& rng) override; double jump(std::mt19937_64 &rng) override;
[[nodiscard]] std::unique_ptr<BaseMotion> clone() const override; [[nodiscard]] std::unique_ptr<BaseMotion> clone() const override;
[[nodiscard]] std::string toString() const override; [[nodiscard]] std::string toString() const override;
}; };
} } // namespace motions
#endif //CONEWOBBLE_H #endif // CONEWOBBLE_H

56
src/motions/coordinates.cpp
View File

@@ -1,41 +1,43 @@
#include "coordinates.h" #include "coordinates.h"
#include <cmath> #include <cmath>
#include <iostream>
namespace coordinates { namespace coordinates {
SphericalPos rotate(const SphericalPos& old_pos, const double alpha, const double beta) { SphericalPos rotate(const SphericalPos &old_pos, const double alpha,
const double sin_alpha{std::sin(alpha)}; const double beta) {
const double cos_alpha{std::cos(alpha)}; const double sin_alpha{std::sin(alpha)};
const double cos_alpha{std::cos(alpha)};
const double sin_beta{std::sin(beta)}; const double sin_beta{std::sin(beta)};
const double cos_beta{std::cos(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) { if (std::abs(cos_theta) == 1) {
return xyz_to_spherical(CartesianPos{cos_alpha * cos_beta, cos_alpha * sin_beta, cos_alpha * cos_theta}); 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{ const auto new_pos = CartesianPos{
cos_alpha * x + sin_alpha * (-cos_beta * y - sin_beta * x * z) / norm, 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 * y + sin_alpha * (cos_beta * x - sin_beta * y * z) / norm,
cos_alpha * z + sin_alpha * sin_beta * 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) { CartesianPos spherical_to_xyz(const SphericalPos &pos) {
const double sin_theta = std::sin(std::acos(pos.cos_theta)); 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}; return {sin_theta * std::cos(pos.phi), sin_theta * std::sin(pos.phi),
} pos.cos_theta};
}
SphericalPos xyz_to_spherical(const CartesianPos& pos) { SphericalPos xyz_to_spherical(const CartesianPos &pos) {
return {pos.z, std::atan2(pos.y, pos.x)}; return {pos.z, std::atan2(pos.y, pos.x)};
} }
} } // namespace coordinates

28
src/motions/coordinates.h
View File

@@ -3,20 +3,20 @@
#define COORDINATES_H #define COORDINATES_H
namespace coordinates { namespace coordinates {
struct SphericalPos { struct SphericalPos {
double cos_theta; double cos_theta;
double phi; double phi;
}; };
struct CartesianPos { struct CartesianPos {
double x; double x;
double y; double y;
double z; double z;
}; };
SphericalPos rotate(const SphericalPos&, double, double); SphericalPos rotate(const SphericalPos &, double, double);
CartesianPos spherical_to_xyz(const SphericalPos&); CartesianPos spherical_to_xyz(const SphericalPos &);
SphericalPos xyz_to_spherical(const CartesianPos&); SphericalPos xyz_to_spherical(const CartesianPos &);
} } // namespace coordinates
#endif //COORDINATES_H #endif // COORDINATES_H

8
src/motions/diffusivemotion.cpp
View File

@@ -1,8 +1,8 @@
#include "diffusivemotion.h" #include "diffusivemotion.h"
namespace motions { namespace motions {
void DiffusiveMotion::initialize(std::mt19937_64& rng) { void DiffusiveMotion::initialize(std::mt19937_64 &rng) {
m_prev_pos = draw_position(rng); m_prev_pos = draw_position(rng);
m_initial_omega = omega_q(m_prev_pos); m_initial_omega = omega_q(m_prev_pos);
}
} }
} // namespace motions

18
src/motions/diffusivemotion.h
View File

@@ -5,15 +5,15 @@
#include "coordinates.h" #include "coordinates.h"
namespace motions { namespace motions {
class DiffusiveMotion : public BaseMotion { class DiffusiveMotion : public BaseMotion {
public: public:
using BaseMotion::BaseMotion; using BaseMotion::BaseMotion;
void initialize(std::mt19937_64& rng) override; void initialize(std::mt19937_64 &rng) override;
protected: protected:
coordinates::SphericalPos m_prev_pos{0., 0.}; coordinates::SphericalPos m_prev_pos{0., 0.};
}; };
} } // namespace motions
#endif //DIFFUSIVEMOTION_H #endif // DIFFUSIVEMOTION_H

62
src/motions/foursitejump.cpp
View File

@@ -1,40 +1,42 @@
#include "foursitejump.h" #include "foursitejump.h"
#include "coordinates.h"
#include "../utils/registry.h" #include "../utils/registry.h"
#include "coordinates.h"
static AutoRegister<motions::BaseMotion, motions::FourSiteTetrahedron> reg("FourSiteTetrahedral"); static AutoRegister<motions::BaseMotion, motions::FourSiteTetrahedron>
reg("FourSiteTetrahedral");
namespace motions { namespace motions {
FourSiteTetrahedron::FourSiteTetrahedron(const double delta, const double eta) : FourSiteTetrahedron::FourSiteTetrahedron(const double delta, const double eta)
BaseMotion(std::string{"FourSiteTetrahedral"}, delta, 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) { void FourSiteTetrahedron::initialize(std::mt19937_64 &rng) {
const auto pos = draw_position(rng); const auto pos = draw_position(rng);
m_corners[0] = omega_q(pos); 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++) { for (int i = 1; i < 4; i++) {
auto corner_pos = rotate(pos, m_beta, alpha + (i-1) * 2*M_PI/3.); auto corner_pos = rotate(pos, m_beta, alpha + (i - 1) * 2 * M_PI / 3.);
m_corners[i] = omega_q(corner_pos); m_corners[i] = omega_q(corner_pos);
} }
m_initial_omega = FourSiteTetrahedron::jump(rng); 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<BaseMotion> FourSiteTetrahedron::clone() const {
return std::make_unique<FourSiteTetrahedron>(*this);
}
std::string FourSiteTetrahedron::toString() const {
return {"FourSiteTetrahedral"};
}
} }
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<BaseMotion> FourSiteTetrahedron::clone() const {
return std::make_unique<FourSiteTetrahedron>(*this);
}
std::string FourSiteTetrahedron::toString() const {
return {"FourSiteTetrahedral"};
}
} // namespace motions

35
src/motions/foursitejump.h
View File

@@ -3,29 +3,28 @@
#include "base.h" #include "base.h"
#include <cmath>
#include <array> #include <array>
#include <cmath>
namespace motions { namespace motions {
class FourSiteTetrahedron final : public BaseMotion { class FourSiteTetrahedron final : public BaseMotion {
public: public:
FourSiteTetrahedron(double, double); FourSiteTetrahedron(double, double);
FourSiteTetrahedron(); FourSiteTetrahedron();
void initialize(std::mt19937_64& rng) override; void initialize(std::mt19937_64 &rng) override;
double jump(std::mt19937_64& rng) override; double jump(std::mt19937_64 &rng) override;
[[nodiscard]] std::unique_ptr<BaseMotion> clone() const override; [[nodiscard]] std::unique_ptr<BaseMotion> clone() const override;
[[nodiscard]] std::string toString() const override; [[nodiscard]] std::string toString() const override;
private: private:
const double m_beta{std::acos(-1/3.)}; const double m_beta{std::acos(-1 / 3.)};
std::array<double, 4> m_corners{}; std::array<double, 4> m_corners{};
int m_corner_idx{0}; int m_corner_idx{0};
std::uniform_int_distribution<> m_chooser{1, 3}; std::uniform_int_distribution<> m_chooser{1, 3};
};
}; } // namespace motions
} #endif // RWSIM_MOTIONTETRAHEDRAL_H
#endif //RWSIM_MOTIONTETRAHEDRAL_H

62
src/motions/isosmallangle.cpp
View File

@@ -1,37 +1,41 @@
#include "isosmallangle.h" #include "isosmallangle.h"
#include "coordinates.h"
#include "../utils/registry.h" #include "../utils/registry.h"
#include "coordinates.h"
static AutoRegister<motions::BaseMotion, motions::SmallAngle> reg("IsotropicAngle"); static AutoRegister<motions::BaseMotion, motions::SmallAngle>
reg("IsotropicAngle");
namespace motions { namespace motions {
SmallAngle::SmallAngle(const double delta, const double eta, const double chi) : SmallAngle::SmallAngle(const double delta, const double eta, const double chi)
DiffusiveMotion(std::string("IsotropicAngle"), delta, eta), m_chi(chi * M_PI / 180.0) {} : DiffusiveMotion(std::string("IsotropicAngle"), delta, eta),
SmallAngle::SmallAngle() : DiffusiveMotion(std::string("IsotropicAngle")) {} m_chi(chi * M_PI / 180.0) {}
SmallAngle::SmallAngle() : DiffusiveMotion(std::string("IsotropicAngle")) {}
double SmallAngle::jump(std::mt19937_64& rng) { double SmallAngle::jump(std::mt19937_64 &rng) {
const double gamma{2 * M_PI * m_uni_dist(rng)}; const double gamma{2 * M_PI * m_uni_dist(rng)};
m_prev_pos = rotate(m_prev_pos, m_chi, gamma); m_prev_pos = rotate(m_prev_pos, m_chi, gamma);
return omega_q(m_prev_pos); return omega_q(m_prev_pos);
}
std::unique_ptr<BaseMotion> SmallAngle::clone() const {
return std::make_unique<SmallAngle>(*this);
}
void SmallAngle::setParameters(const std::unordered_map<std::string, double> &parameters) {
m_chi = parameters.at("angle") * M_PI / 180.0;
BaseMotion::setParameters(parameters);
}
std::unordered_map<std::string, double> 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)};
}
} }
std::unique_ptr<BaseMotion> SmallAngle::clone() const {
return std::make_unique<SmallAngle>(*this);
}
void SmallAngle::setParameters(
const std::unordered_map<std::string, double> &parameters) {
m_chi = parameters.at("angle") * M_PI / 180.0;
BaseMotion::setParameters(parameters);
}
std::unordered_map<std::string, double> 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

29
src/motions/isosmallangle.h
View File

@@ -5,19 +5,20 @@
#include "diffusivemotion.h" #include "diffusivemotion.h"
namespace motions { namespace motions {
class SmallAngle final : public DiffusiveMotion { class SmallAngle final : public DiffusiveMotion {
public: public:
SmallAngle(double, double, double); SmallAngle(double, double, double);
SmallAngle(); SmallAngle();
double jump(std::mt19937_64& rng) override; double jump(std::mt19937_64 &rng) override;
[[nodiscard]] std::unique_ptr<BaseMotion> clone() const override; [[nodiscard]] std::unique_ptr<BaseMotion> clone() const override;
void setParameters(const std::unordered_map<std::string, double> &) override; void setParameters(const std::unordered_map<std::string, double> &) override;
[[nodiscard]] std::unordered_map<std::string, double> getParameters() const override; [[nodiscard]] std::unordered_map<std::string, double>
[[nodiscard]] std::string toString() const override; getParameters() const override;
[[nodiscard]] std::string toString() const override;
private: private:
double m_chi{0}; double m_chi{0};
}; };
} } // namespace motions
#endif //RWSIM_MOTIONISOSMALLANGLE_H #endif // RWSIM_MOTIONISOSMALLANGLE_H

84
src/motions/nsiteconejump.cpp
View File

@@ -1,61 +1,63 @@
#include "nsiteconejump.h" #include "nsiteconejump.h"
#include "coordinates.h"
#include "../utils/registry.h" #include "../utils/registry.h"
#include "coordinates.h"
#include <cmath> #include <cmath>
#include <vector> #include <vector>
static AutoRegister<motions::BaseMotion, motions::NSiteJumpOnCone> reg("NSiteConeJump"); static AutoRegister<motions::BaseMotion, motions::NSiteJumpOnCone>
reg("NSiteConeJump");
namespace motions { namespace motions {
NSiteJumpOnCone::NSiteJumpOnCone(const double delta, const double eta, const int num_sites, const double chi) : NSiteJumpOnCone::NSiteJumpOnCone(const double delta, const double eta,
BaseMotion("NSiteJumpOnCone", delta, eta), const int num_sites, const double chi)
m_num_sites(num_sites), : BaseMotion("NSiteJumpOnCone", delta, eta), m_num_sites(num_sites),
m_chi(chi*M_PI/180.) {} m_chi(chi * M_PI / 180.) {}
NSiteJumpOnCone::NSiteJumpOnCone() : BaseMotion("NSiteJumpOnCone") { } NSiteJumpOnCone::NSiteJumpOnCone() : BaseMotion("NSiteJumpOnCone") {}
void NSiteJumpOnCone::initialize(std::mt19937_64& rng) { void NSiteJumpOnCone::initialize(std::mt19937_64 &rng) {
m_sites = std::vector<double>(m_num_sites); m_sites = std::vector<double>(m_num_sites);
m_chooser = std::uniform_int_distribution<>{1, m_num_sites - 1}; 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 alpha = m_uni_dist(rng) * 2 * M_PI;
const double steps = 2*M_PI / m_num_sites; const double steps = 2 * M_PI / m_num_sites;
for (int i = 0; i < m_num_sites; i++) { for (int i = 0; i < m_num_sites; i++) {
m_sites[i] = omega_q(rotate(m_axis, m_chi, i * steps + alpha)); m_sites[i] = omega_q(rotate(m_axis, m_chi, i * steps + alpha));
} }
} }
double NSiteJumpOnCone::jump(std::mt19937_64& rng) { double NSiteJumpOnCone::jump(std::mt19937_64 &rng) {
m_cone_idx += m_chooser(rng); m_cone_idx += m_chooser(rng);
m_cone_idx %= m_num_sites; m_cone_idx %= m_num_sites;
return m_sites[m_cone_idx]; return m_sites[m_cone_idx];
} }
std::unique_ptr<BaseMotion> NSiteJumpOnCone::clone() const { std::unique_ptr<BaseMotion> NSiteJumpOnCone::clone() const {
return std::make_unique<NSiteJumpOnCone>(*this); return std::make_unique<NSiteJumpOnCone>(*this);
} }
void NSiteJumpOnCone::setParameters(const std::unordered_map<std::string, double> &parameters) { void NSiteJumpOnCone::setParameters(
BaseMotion::setParameters(parameters); const std::unordered_map<std::string, double> &parameters) {
m_num_sites = static_cast<int>(parameters.at("num_sites")); BaseMotion::setParameters(parameters);
m_chi = parameters.at("chi") * M_PI/180.; m_num_sites = static_cast<int>(parameters.at("num_sites"));
} m_chi = parameters.at("chi") * M_PI / 180.;
}
std::unordered_map<std::string, double> NSiteJumpOnCone::getParameters() const { std::unordered_map<std::string, double> NSiteJumpOnCone::getParameters() const {
auto parameter = BaseMotion::getParameters(); auto parameter = BaseMotion::getParameters();
parameter["num_sites"] = m_num_sites; parameter["num_sites"] = m_num_sites;
parameter["chi"] = m_chi * 180. / M_PI; parameter["chi"] = m_chi * 180. / M_PI;
return parameter; return parameter;
} }
std::string NSiteJumpOnCone::toString() const { std::string NSiteJumpOnCone::toString() const {
return std::to_string(m_num_sites) + "SiteJumpOnCone/angle=" + std::to_string(m_chi*180/M_PI); return std::to_string(m_num_sites) +
} "SiteJumpOnCone/angle=" + std::to_string(m_chi * 180 / M_PI);
}
} // motions } // namespace motions

41
src/motions/nsiteconejump.h
View File

@@ -6,28 +6,29 @@
#include <vector> #include <vector>
namespace motions { namespace motions {
class NSiteJumpOnCone final : public BaseMotion { class NSiteJumpOnCone final : public BaseMotion {
public: public:
NSiteJumpOnCone(double, double, int, double); NSiteJumpOnCone(double, double, int, double);
NSiteJumpOnCone(); NSiteJumpOnCone();
void initialize(std::mt19937_64& rng) override; void initialize(std::mt19937_64 &rng) override;
double jump(std::mt19937_64& rng) override; double jump(std::mt19937_64 &rng) override;
[[nodiscard]] std::unique_ptr<BaseMotion> clone() const override; [[nodiscard]] std::unique_ptr<BaseMotion> clone() const override;
[[nodiscard]] std::string toString() const override; [[nodiscard]] std::string toString() const override;
void setParameters(const std::unordered_map<std::string, double> &) override; void setParameters(const std::unordered_map<std::string, double> &) override;
[[nodiscard]] std::unordered_map<std::string, double> getParameters() const override; [[nodiscard]] std::unordered_map<std::string, double>
getParameters() const override;
private: private:
int m_num_sites{2}; int m_num_sites{2};
std::vector<double> m_sites{}; std::vector<double> m_sites{};
double m_chi{1./2.}; double m_chi{1. / 2.};
int m_cone_idx = 0; int m_cone_idx = 0;
coordinates::SphericalPos m_axis{1, 0}; coordinates::SphericalPos m_axis{1, 0};
std::uniform_int_distribution<> m_chooser; std::uniform_int_distribution<> m_chooser;
}; };
} // motions } // namespace motions
#endif //NSITEJUMPONCONE_H #endif // NSITEJUMPONCONE_H

32
src/motions/random.cpp
View File

@@ -1,27 +1,25 @@
#include "random.h" #include "random.h"
#include "../utils/registry.h" #include "../utils/registry.h"
static AutoRegister<motions::BaseMotion, motions::RandomJump> reg("RandomJump"); static AutoRegister<motions::BaseMotion, motions::RandomJump> reg("RandomJump");
namespace motions { 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 { std::string RandomJump::toString() const { return {"RandomJump"}; }
return {"RandomJump"};
}
std::unique_ptr<BaseMotion> RandomJump::clone() const { std::unique_ptr<BaseMotion> RandomJump::clone() const {
return std::make_unique<RandomJump>(*this); return std::make_unique<RandomJump>(*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));
}
} }
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

23
src/motions/random.h
View File

@@ -1,21 +1,20 @@
#ifndef RWSIM_MOTIONRANDOMJUMP_H #ifndef RWSIM_MOTIONRANDOMJUMP_H
#define RWSIM_MOTIONRANDOMJUMP_H #define RWSIM_MOTIONRANDOMJUMP_H
#include "base.h" #include "base.h"
namespace motions { namespace motions {
class RandomJump final : public BaseMotion { class RandomJump final : public BaseMotion {
public: public:
RandomJump(double, double); RandomJump(double, double);
RandomJump(); RandomJump();
[[nodiscard]] std::string toString() const override; [[nodiscard]] std::string toString() const override;
[[nodiscard]] std::unique_ptr<BaseMotion> clone() const override; [[nodiscard]] std::unique_ptr<BaseMotion> clone() const override;
void initialize(std::mt19937_64& rng) override; void initialize(std::mt19937_64 &rng) override;
double jump(std::mt19937_64& rng) override; double jump(std::mt19937_64 &rng) override;
}; };
} } // namespace motions
#endif //RWSIM_MOTIONRANDOMJUMP_H #endif // RWSIM_MOTIONRANDOMJUMP_H

40
src/motions/rjoac.cpp
View File

@@ -1,24 +1,28 @@
#include "rjoac.h" #include "rjoac.h"
#include "coordinates.h"
#include "../utils/registry.h" #include "../utils/registry.h"
#include "coordinates.h"
static AutoRegister<motions::BaseMotion, motions::RandomJumpOnCone> reg("RandomJumpOnCone"); static AutoRegister<motions::BaseMotion, motions::RandomJumpOnCone>
reg("RandomJumpOnCone");
namespace motions { 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(const double delta, const double eta,
RandomJumpOnCone::RandomJumpOnCone() : ConeMotion("RJ on a Cone") {} const double chi)
: ConeMotion("RJ on a Cone", delta, eta) {
double RandomJumpOnCone::jump(std::mt19937_64& rng) { m_angle = chi;
const double phi = 2 * M_PI * m_uni_dist(rng);
return omega_q(rotate(m_axis, m_angle, phi));
}
std::unique_ptr<BaseMotion> RandomJumpOnCone::clone() const {
return std::make_unique<RandomJumpOnCone>(*this);
}
std::string RandomJumpOnCone::toString() const {
return std::string("RandomJumpOnCone/angle=") + std::to_string(m_angle);
}
} }
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<BaseMotion> RandomJumpOnCone::clone() const {
return std::make_unique<RandomJumpOnCone>(*this);
}
std::string RandomJumpOnCone::toString() const {
return std::string("RandomJumpOnCone/angle=") + std::to_string(m_angle);
}
} // namespace motions

20
src/motions/rjoac.h
View File

@@ -4,15 +4,15 @@
#include "conemotion.h" #include "conemotion.h"
namespace motions { namespace motions {
class RandomJumpOnCone final: public ConeMotion { class RandomJumpOnCone final : public ConeMotion {
public: public:
RandomJumpOnCone(double, double, double); RandomJumpOnCone(double, double, double);
RandomJumpOnCone(); RandomJumpOnCone();
double jump(std::mt19937_64& rng) override; double jump(std::mt19937_64 &rng) override;
[[nodiscard]] std::unique_ptr<BaseMotion> clone() const override; [[nodiscard]] std::unique_ptr<BaseMotion> clone() const override;
[[nodiscard]] std::string toString() const override; [[nodiscard]] std::string toString() const override;
}; };
} } // namespace motions
#endif //RJOAC_H #endif // RJOAC_H

68
src/motions/sixsitejump.cpp
View File

@@ -1,43 +1,47 @@
#include "sixsitejump.h" #include "sixsitejump.h"
#include "coordinates.h"
#include "../utils/registry.h" #include "../utils/registry.h"
#include "coordinates.h"
static AutoRegister<motions::BaseMotion, motions::SixSiteOctahedronC3> reg("SixSiteOctahedralC3"); static AutoRegister<motions::BaseMotion, motions::SixSiteOctahedronC3>
reg("SixSiteOctahedralC3");
namespace motions { namespace motions {
SixSiteOctahedronC3::SixSiteOctahedronC3(const double delta, const double eta, const double chi) : SixSiteOctahedronC3::SixSiteOctahedronC3(const double delta, const double eta,
BaseMotion(std::string{"SixSiteOctahedral"}, delta, eta), const double chi)
m_chi{chi*M_PI/180.} {} : 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) { void SixSiteOctahedronC3::initialize(std::mt19937_64 &rng) {
const coordinates::SphericalPos c3_axis = draw_position(rng); const coordinates::SphericalPos c3_axis = draw_position(rng);
const double alpha = 2. * M_PI * m_uni_dist(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++) { 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] =
m_corners[2*i+1] = omega_q(rotate(c3_axis, m_chi_opposite, alpha + i * 2./3.*M_PI + M_PI/3.)); 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); 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<BaseMotion> SixSiteOctahedronC3::clone() const {
return std::make_unique<SixSiteOctahedronC3>(*this);
}
std::string SixSiteOctahedronC3::toString() const {
return {"SixSiteOctahedral/angle=" + std::to_string(m_chi / M_PI * 180.)};
}
} }
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<BaseMotion> SixSiteOctahedronC3::clone() const {
return std::make_unique<SixSiteOctahedronC3>(*this);
}
std::string SixSiteOctahedronC3::toString() const {
return {"SixSiteOctahedral/angle=" + std::to_string(m_chi / M_PI * 180.)};
}
} // namespace motions

36
src/motions/sixsitejump.h
View File

@@ -1,32 +1,30 @@
#ifndef SIXSITEJUMP_H #ifndef SIXSITEJUMP_H
#define SIXSITEJUMP_H #define SIXSITEJUMP_H
#include "base.h" #include "base.h"
#include <cmath>
#include <array> #include <array>
#include <cmath>
namespace motions { namespace motions {
class SixSiteOctahedronC3 final : public BaseMotion { class SixSiteOctahedronC3 final : public BaseMotion {
public: public:
SixSiteOctahedronC3(double, double, double); SixSiteOctahedronC3(double, double, double);
SixSiteOctahedronC3(); SixSiteOctahedronC3();
void initialize(std::mt19937_64& rng) override; void initialize(std::mt19937_64 &rng) override;
double jump(std::mt19937_64& rng) override; double jump(std::mt19937_64 &rng) override;
[[nodiscard]] std::unique_ptr<BaseMotion> clone() const override; [[nodiscard]] std::unique_ptr<BaseMotion> clone() const override;
[[nodiscard]] std::string toString() const override; [[nodiscard]] std::string toString() const override;
private: private:
double m_chi{std::acos(-1.0 / 3.0)}; // 54.74 deg double m_chi{std::acos(-1.0 / 3.0)}; // 54.74 deg
std::array<double, 6> m_corners{}; std::array<double, 6> m_corners{};
int m_corner_idx{0}; int m_corner_idx{0};
std::uniform_int_distribution<> m_chooser{1, 5}; std::uniform_int_distribution<> m_chooser{1, 5};
};
}; } // namespace motions
} #endif // SIXSITEJUMP_H
#endif //SIXSITEJUMP_H

189
src/sims.cpp
View File

@@ -2,131 +2,124 @@
#include "utils/functions.h" #include "utils/functions.h"
#include "utils/io.h" #include "utils/io.h"
#include <iostream>
#include <chrono> #include <chrono>
#include <iostream>
#include <omp.h> #include <omp.h>
void run_simulation(Experiment &experiment,
const std::unordered_map<std::string, double> &parameter,
const std::unordered_map<std::string, double> &optional,
Dynamics &dynamics, std::mt19937_64 &rng) {
const int num_walker = static_cast<int>(parameter.at("num_walker"));
void run_simulation( dynamics.setParameters(parameter);
Experiment& experiment, experiment.setup(parameter, optional);
const std::unordered_map<std::string, double>& parameter,
const std::unordered_map<std::string, double>& optional,
Dynamics& dynamics,
std::mt19937_64& rng
) {
const int num_walker = static_cast<int>(parameter.at("num_walker"));
dynamics.setParameters(parameter); const auto start = printStart(optional);
experiment.setup(parameter, 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<std::mt19937_64> 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 // Create per-thread clones of dynamics and experiment
std::vector<std::mt19937_64> thread_rngs; std::vector<std::unique_ptr<Dynamics>> thread_dynamics;
thread_rngs.reserve(num_threads); std::vector<std::unique_ptr<Experiment>> thread_experiments;
for (int i = 0; i < num_threads; i++) { for (int i = 0; i < num_threads; i++) {
thread_rngs.emplace_back(rng()); 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 // Merge per-thread results
std::vector<std::unique_ptr<Dynamics>> thread_dynamics; for (int i = 0; i < num_threads; i++) {
std::vector<std::unique_ptr<Experiment>> thread_experiments; experiment.merge(*thread_experiments[i]);
for (int i = 0; i < num_threads; i++) { }
thread_dynamics.push_back(dynamics.clone());
thread_experiments.push_back(experiment.clone());
}
int steps_done = 0; const auto directory = make_directory(dynamics.toString());
auto last_print_out = std::chrono::system_clock::now(); experiment.save(directory);
printEnd(start);
#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);
} }
Trajectory make_trajectory(Dynamics &dynamics, const double t_max,
std::mt19937_64 &rng) {
double t_passed = 0;
double phase = 0;
Trajectory make_trajectory( dynamics.initialize(rng);
Dynamics& dynamics,
const double t_max,
std::mt19937_64& rng
) {
double t_passed = 0;
double phase = 0;
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.time.emplace_back(t_passed);
traj.phase.emplace_back(phase); traj.phase.emplace_back(phase);
traj.omega.emplace_back(omega); traj.omega.emplace_back(omega);
}
while (t_passed < t_max) { return traj;
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;
} }
std::chrono::system_clock::time_point
printStart(const std::unordered_map<std::string, double> &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<std::string, double> &optional) { std::cout << "Random walk for ";
const auto start = std::chrono::system_clock::now(); for (const auto &[key, value] : optional) {
const time_t start_time = std::chrono::system_clock::to_time_t(start); std::cout << key << " = " << value << "; ";
}
std::cout << std::endl;
std::cout << "Start: " << ctime(&start_time);
std::cout << "Random walk for "; return start;
for (const auto& [key, value] : optional) {
std::cout << key << " = " << value << "; ";
}
std::cout << std::endl;
std::cout << "Start: " << ctime(&start_time);
return start;
} }
void printEnd(const std::chrono::system_clock::time_point 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<float> duration = end - start; const std::chrono::duration<float> duration = end - start;
const time_t end_time = std::chrono::system_clock::to_time_t(end); const time_t end_time = std::chrono::system_clock::to_time_t(end);
std::cout << "End: " << ctime(&end_time); std::cout << "End: " << ctime(&end_time);
std::cout << "Duration: " << duration.count() << "s\n" << std::endl; std::cout << "Duration: " << duration.count() << "s\n" << std::endl;
} }

22
src/sims.h
View File

@@ -4,21 +4,21 @@
#include "dynamics/base.h" #include "dynamics/base.h"
#include "experiments/base.h" #include "experiments/base.h"
#include <unordered_map>
#include <string>
#include <chrono> #include <chrono>
#include <random> #include <random>
#include <string>
#include <unordered_map>
void run_simulation( void run_simulation(Experiment &experiment,
Experiment& experiment, const std::unordered_map<std::string, double> &parameter,
const std::unordered_map<std::string, double>& parameter, const std::unordered_map<std::string, double> &optional,
const std::unordered_map<std::string, double>& optional, Dynamics &dynamics, std::mt19937_64 &rng);
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<std::string, double> &optional); std::chrono::system_clock::time_point
printStart(const std::unordered_map<std::string, double> &optional);
void printEnd(std::chrono::system_clock::time_point start); void printEnd(std::chrono::system_clock::time_point start);
#endif //RWSIM_SIMS_H #endif // RWSIM_SIMS_H

12
src/times/CMakeLists.txt
View File

@@ -1,9 +1,5 @@
add_library( add_library(RWTime OBJECT
RWTime OBJECT base.cpp base.h
base.cpp delta.cpp delta.h
base.h lognormal.cpp lognormal.h
delta.cpp
delta.h
lognormal.cpp
lognormal.h
) )

44
src/times/base.cpp
View File

@@ -2,26 +2,30 @@
#include "../utils/registry.h" #include "../utils/registry.h"
namespace times { 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 { double BaseDistribution::tau_wait(std::mt19937_64 &rng) const {
return std::exponential_distribution(1./m_tau_jump)(rng); return std::exponential_distribution(1. / m_tau_jump)(rng);
}
void BaseDistribution::setParameters(const std::unordered_map<std::string, double> &parameters) {
m_tau = parameters.at("tau");
}
std::unordered_map<std::string, double> BaseDistribution::getParameters() const {
return std::unordered_map<std::string, double>{
{"tau", m_tau},
};
}
std::unique_ptr<BaseDistribution> BaseDistribution::createFromInput(const std::string& input) {
return Registry<BaseDistribution>::create(input);
}
} }
void BaseDistribution::setParameters(
const std::unordered_map<std::string, double> &parameters) {
m_tau = parameters.at("tau");
}
std::unordered_map<std::string, double>
BaseDistribution::getParameters() const {
return std::unordered_map<std::string, double>{
{"tau", m_tau},
};
}
std::unique_ptr<BaseDistribution>
BaseDistribution::createFromInput(const std::string &input) {
return Registry<BaseDistribution>::create(input);
}
} // namespace times

46
src/times/base.h
View File

@@ -6,33 +6,35 @@
#include <unordered_map> #include <unordered_map>
namespace times { namespace times {
class BaseDistribution { class BaseDistribution {
public: public:
virtual ~BaseDistribution() = default; virtual ~BaseDistribution() = default;
BaseDistribution(std::string, double); BaseDistribution(std::string, double);
explicit BaseDistribution(std::string); explicit BaseDistribution(std::string);
[[nodiscard]] double getTau() const { return m_tau; } [[nodiscard]] double getTau() const { return m_tau; }
void setTau(const double tau) { m_tau = tau; } void setTau(const double tau) { m_tau = tau; }
[[nodiscard]] std::string getName() const { return m_name; } [[nodiscard]] std::string getName() const { return m_name; }
virtual void setParameters(const std::unordered_map<std::string, double>&); virtual void setParameters(const std::unordered_map<std::string, double> &);
[[nodiscard]] virtual std::unordered_map<std::string, double> getParameters() const; [[nodiscard]] virtual std::unordered_map<std::string, double>
getParameters() const;
virtual void initialize(std::mt19937_64& rng) = 0; virtual void initialize(std::mt19937_64 &rng) = 0;
[[nodiscard]] double tau_wait(std::mt19937_64& rng) const; [[nodiscard]] double tau_wait(std::mt19937_64 &rng) const;
[[nodiscard]] virtual std::unique_ptr<BaseDistribution> clone() const = 0; [[nodiscard]] virtual std::unique_ptr<BaseDistribution> clone() const = 0;
[[nodiscard]] virtual std::string toString() const = 0; [[nodiscard]] virtual std::string toString() const = 0;
static std::unique_ptr<BaseDistribution> createFromInput(const std::string& input); static std::unique_ptr<BaseDistribution>
createFromInput(const std::string &input);
protected: protected:
std::string m_name{"BaseDistribution"}; std::string m_name{"BaseDistribution"};
double m_tau{1.}; double m_tau{1.};
double m_tau_jump{1.}; double m_tau_jump{1.};
}; };
} } // namespace times
#endif //RWSIM_TIMESBASE_H #endif // RWSIM_TIMESBASE_H

27
src/times/delta.cpp
View File

@@ -1,21 +1,22 @@
#include "delta.h" #include "delta.h"
#include "../utils/registry.h" #include "../utils/registry.h"
static AutoRegister<times::BaseDistribution, times::DeltaDistribution> reg("Delta"); static AutoRegister<times::BaseDistribution, times::DeltaDistribution>
reg("Delta");
namespace times { namespace times {
DeltaDistribution::DeltaDistribution(const double tau) : BaseDistribution(std::string("Delta"), tau) {} DeltaDistribution::DeltaDistribution(const double tau)
DeltaDistribution::DeltaDistribution() : BaseDistribution(std::string("Delta")) {} : BaseDistribution(std::string("Delta"), tau) {}
DeltaDistribution::DeltaDistribution()
: BaseDistribution(std::string("Delta")) {}
void DeltaDistribution::initialize(std::mt19937_64&) { void DeltaDistribution::initialize(std::mt19937_64 &) { m_tau_jump = m_tau; }
m_tau_jump = m_tau;
}
std::unique_ptr<BaseDistribution> DeltaDistribution::clone() const { std::unique_ptr<BaseDistribution> DeltaDistribution::clone() const {
return std::make_unique<DeltaDistribution>(*this); return std::make_unique<DeltaDistribution>(*this);
}
std::string DeltaDistribution::toString() const {
return {"Delta/tau=" + std::to_string(m_tau)};
}
} }
std::string DeltaDistribution::toString() const {
return {"Delta/tau=" + std::to_string(m_tau)};
}
} // namespace times

20
src/times/delta.h
View File

@@ -4,15 +4,15 @@
#include "base.h" #include "base.h"
namespace times { namespace times {
class DeltaDistribution final : public BaseDistribution { class DeltaDistribution final : public BaseDistribution {
public: public:
explicit DeltaDistribution(double); explicit DeltaDistribution(double);
DeltaDistribution(); DeltaDistribution();
void initialize(std::mt19937_64& rng) override; void initialize(std::mt19937_64 &rng) override;
[[nodiscard]] std::unique_ptr<BaseDistribution> clone() const override; [[nodiscard]] std::unique_ptr<BaseDistribution> clone() const override;
[[nodiscard]] std::string toString() const override; [[nodiscard]] std::string toString() const override;
}; };
} } // namespace times
#endif //RWSIM_TIMESDELTA_H #endif // RWSIM_TIMESDELTA_H

60
src/times/lognormal.cpp
View File

@@ -2,33 +2,41 @@
#include "../utils/registry.h" #include "../utils/registry.h"
#include <cmath> #include <cmath>
static AutoRegister<times::BaseDistribution, times::LogNormalDistribution> reg("LogNormal"); static AutoRegister<times::BaseDistribution, times::LogNormalDistribution>
reg("LogNormal");
namespace times { 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(const double tau,
LogNormalDistribution::LogNormalDistribution() : BaseDistribution(std::string("LogNormal")) {} 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<std::string, double> &parameters) { void LogNormalDistribution::setParameters(
m_sigma = parameters.at("sigma"); const std::unordered_map<std::string, double> &parameters) {
BaseDistribution::setParameters(parameters); m_sigma = parameters.at("sigma");
} BaseDistribution::setParameters(parameters);
std::unordered_map<std::string, double> 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<BaseDistribution> LogNormalDistribution::clone() const {
return std::make_unique<LogNormalDistribution>(*this);
}
std::string LogNormalDistribution::toString() const {
return {"LogNormal/tau=" + std::to_string(m_tau) + "/sigma=" + std::to_string(m_sigma)};
}
} }
std::unordered_map<std::string, double>
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<BaseDistribution> LogNormalDistribution::clone() const {
return std::make_unique<LogNormalDistribution>(*this);
}
std::string LogNormalDistribution::toString() const {
return {"LogNormal/tau=" + std::to_string(m_tau) +
"/sigma=" + std::to_string(m_sigma)};
}
} // namespace times

31
src/times/lognormal.h
View File

@@ -5,23 +5,24 @@
#include <random> #include <random>
namespace times { namespace times {
class LogNormalDistribution final : public BaseDistribution { class LogNormalDistribution final : public BaseDistribution {
public: public:
LogNormalDistribution(double, double); LogNormalDistribution(double, double);
LogNormalDistribution(); LogNormalDistribution();
void setParameters(const std::unordered_map<std::string, double> &) override; void setParameters(const std::unordered_map<std::string, double> &) override;
[[nodiscard]] std::unordered_map<std::string, double> getParameters() const override; [[nodiscard]] std::unordered_map<std::string, double>
getParameters() const override;
[[nodiscard]] std::string toString() const override; [[nodiscard]] std::string toString() const override;
[[nodiscard]] std::unique_ptr<BaseDistribution> clone() const override; [[nodiscard]] std::unique_ptr<BaseDistribution> clone() const override;
void initialize(std::mt19937_64& rng) override; void initialize(std::mt19937_64 &rng) override;
private: private:
double m_sigma{1}; double m_sigma{1};
std::lognormal_distribution<> m_distribution; std::lognormal_distribution<> m_distribution;
}; };
} } // namespace times
#endif //LOGNORMAL_H #endif // LOGNORMAL_H

71
src/utils/functions.cpp
View File

@@ -1,51 +1,54 @@
#include <vector>
#include <chrono> #include <chrono>
#include <iostream> #include <iostream>
#include <vector>
int nearest_index(const std::vector<double> &x_ref, const double x,
int nearest_index(const std::vector<double> &x_ref, const double x, int start=0) { int start = 0) {
const int last = static_cast<int>(x_ref.size()) - 2; const int last = static_cast<int>(x_ref.size()) - 2;
while (start < last && x > x_ref[start+1]) { while (start < last && x > x_ref[start + 1]) {
start++; start++;
} }
return start; return start;
} }
double lerp(const std::vector<double>& x_ref, const std::vector<double>& y_ref, const double x, const int i) { double lerp(const std::vector<double> &x_ref, const std::vector<double> &y_ref,
/* const double x, const int i) {
* Linear interpolation between two /*
*/ * Linear interpolation between two
const double x_left = x_ref[i]; */
const double y_left = y_ref[i]; const double x_left = x_ref[i];
const double x_right = x_ref[i+1]; const double y_left = y_ref[i];
const double y_right = y_ref[i+1]; 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<std::chrono::system_clock> printSteps( std::chrono::time_point<std::chrono::system_clock> 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<std::chrono::system_clock> last_print_out, const std::chrono::time_point<std::chrono::system_clock>
last_print_out,
const std::chrono::time_point<std::chrono::system_clock> start, const std::chrono::time_point<std::chrono::system_clock> start,
const int total, const int total, const int steps) {
const int steps const auto now = std::chrono::high_resolution_clock::now();
) {
const auto now = std::chrono::high_resolution_clock::now();
if (const std::chrono::duration<float> duration = now - last_print_out; duration.count() < 10.) { if (const std::chrono::duration<float> duration = now - last_print_out;
return last_print_out; duration.count() < 10.) {
} return last_print_out;
}
const std::chrono::duration<float> duration = now - start; const std::chrono::duration<float> duration = now - start;
const auto passed = duration.count(); const auto passed = duration.count();
std::cout << steps << " of " << total << " steps: " << passed << "s passed; ~"
<< passed * static_cast<float>(total - steps) /
static_cast<float>(steps)
<< "s remaining\n";
std::cout << steps << " of " << total << " steps: " << passed << "s passed; ~" << passed * static_cast<float>(total-steps) / static_cast<float>(steps) << "s remaining\n"; return now;
}
return now;
}

12
src/utils/functions.h
View File

@@ -1,14 +1,16 @@
#ifndef RWSIM_FUNCTIONS_H #ifndef RWSIM_FUNCTIONS_H
#define RWSIM_FUNCTIONS_H #define RWSIM_FUNCTIONS_H
#include <vector>
#include <chrono> #include <chrono>
#include <vector>
int nearest_index(const std::vector<double> &, double, int);
int nearest_index(const std::vector<double>&, double, int); double lerp(const std::vector<double> &, const std::vector<double> &, double,
int);
double lerp(const std::vector<double>&, const std::vector<double>&, double, int); std::chrono::time_point<std::chrono::system_clock>
printSteps(std::chrono::time_point<std::chrono::system_clock>,
std::chrono::time_point<std::chrono::system_clock> printSteps(std::chrono::time_point<std::chrono::system_clock>, std::chrono::time_point<std::chrono::system_clock>, int, int); std::chrono::time_point<std::chrono::system_clock>, int, int);
#endif #endif

368
src/utils/io.cpp
View File

@@ -1,28 +1,28 @@
#include "io.h" #include "io.h"
#include <sstream>
#include <fstream>
#include <iostream>
#include <algorithm> #include <algorithm>
#include <complex>
#include <vector>
#include <iomanip>
#include <map>
#include <string>
#include <filesystem> #include <filesystem>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <map>
#include <sstream>
#include <string>
#include <vector>
std::pair<std::string, double>
get_optional_parameter(std::vector<std::string>::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<std::string, double> get_optional_parameter(std::vector<std::string>::const_iterator &it) { return std::make_pair(stripped_arg, stripped_value);
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);
} }
/** /**
@@ -31,201 +31,197 @@ std::pair<std::string, double> get_optional_parameter(std::vector<std::string>::
* @param argv List of command-line arguments * @param argv List of command-line arguments
* @return Arguments * @return Arguments
*/ */
Arguments parse_args(const int argc, char* argv[]) { Arguments parse_args(const int argc, char *argv[]) {
if (argc < 3) { if (argc < 3) {
throw std::runtime_error("Not enough arguments: missing parameter file"); throw std::runtime_error("Not enough arguments: missing parameter file");
}
Arguments args;
// convert to vector to use iterator for loop
const std::vector<std::string> 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; // Two positional parameters are defined: 1. Location of config file; 2.
// Name of motion model
// convert to vector to use iterator for loop if (args.parameter_file.empty()) {
const std::vector<std::string> input_args(argv + 1, argv + argc); args.parameter_file = *it;
continue;
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");
} }
if (args.parameter_file.empty() || args.motion_type.empty() || args.distribution_type.empty()) { if (args.motion_type.empty()) {
throw std::invalid_argument("Missing argument"); 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<std::string, double>
read_parameter(const std::filesystem::path &infile) {
if (!exists(infile)) {
std::cerr << "File " << infile << " does not exist" << std::endl;
exit(1);
}
std::unordered_map<std::string, double> read_parameter(const std::filesystem::path& infile) { std::ifstream instream(infile);
if (!exists(infile)) {
std::cerr << "File " << infile << " does not exist" << std::endl;
exit(1);
}
std::ifstream instream(infile); std::unordered_map<std::string, double> parameter;
std::unordered_map<std::string, double> parameter; std::string line;
std::string delim = "=";
std::string key;
std::string value;
size_t delim_pos;
std::string line; while (std::getline(instream, line)) {
std::string delim = "="; // skip comment lines starting with #, and empty lines
std::string key; if (line[0] == '#' || line.length() == 1)
std::string value; continue;
size_t delim_pos;
while (std::getline(instream, line)) { // strip spaces from line to have always key=value
// skip comment lines starting with #, and empty lines line.erase(std::remove(line.begin(), line.end(), ' '), line.end());
if (line[0] == '#' || line.length() == 1) continue;
// strip spaces from line to have always key=value // split at '=' character and add to map
line.erase(std::remove(line.begin(), line.end(), ' '), line.end()); 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 return parameter;
delim_pos = line.find('=');
key = line.substr(0, delim_pos);
value = line.substr(delim_pos+1);
parameter[key] = std::stod(value);
}
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) { return path_name;
if (!std::filesystem::create_directories(path_name)) {
std::cout << "Created directory " << path_name << std::endl;
}
return path_name;
} }
void save_parameter_to_file( void save_parameter_to_file(
const std::string& resulttype, const std::string &resulttype, const std::string &directory,
const std::string& directory, const std::unordered_map<std::string, double> &parameter,
const std::unordered_map<std::string, double>& parameter, const std::unordered_map<std::string, double> &optional) {
const std::unordered_map<std::string, double>& optional std::ostringstream parameter_filename;
) { parameter_filename << directory << "/" << resulttype;
std::ostringstream parameter_filename;
parameter_filename << directory << "/" << resulttype;
parameter_filename << std::setprecision(6) << std::scientific; parameter_filename << std::setprecision(6) << std::scientific;
for (const auto& [key, value]: optional) { for (const auto &[key, value] : optional) {
parameter_filename << "_" << key << "=" << value; parameter_filename << "_" << key << "=" << value;
} }
parameter_filename << "_parameter.txt"; parameter_filename << "_parameter.txt";
{ {
// write data to file, columns are secondary axis (echo delay, evolution times) // write data to file, columns are secondary axis (echo delay, evolution
std::string datafile = parameter_filename.str(); // times)
std::ofstream filestream(datafile, std::ios::out); std::string datafile = parameter_filename.str();
std::ofstream filestream(datafile, std::ios::out);
for (const auto& [key, value]: parameter) {
filestream << key << "=" << value << "\n"; 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<double>& x,
const std::map<double, std::vector<double>>& y,
const std::unordered_map<std::string, double>& 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";
}
} }
}
} }
void save_data_to_file( void save_data_to_file(
const std::string& resulttype, const std::string &resulttype, const std::string &directory,
const std::string& directory, const std::vector<double> &x,
const std::vector<double>& x, const std::map<double, std::vector<double>> &y,
const std::vector<double>& y, const std::unordered_map<std::string, double> &optional) {
const std::unordered_map<std::string, double>& optional // make file name
) { std::ostringstream datafile_name;
// make file name datafile_name << directory << "/" << resulttype;
std::ostringstream datafile_name;
datafile_name << directory << "/" << resulttype;
datafile_name << std::setprecision(6) << std::scientific; datafile_name << std::setprecision(6) << std::scientific;
for (const auto& [key, value]: optional) { for (const auto &[key, value] : optional) {
datafile_name << "_" << key << "=" << value; 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 values to file
// write data to file, columns are secondary axis (echo delay, evolution times) auto size = x.size();
std::string datafile = datafile_name.str(); for (unsigned int i = 0; i < size; i++) {
std::ofstream filestream(datafile, std::ios::out); filestream << x[i];
for (const auto &[_, fid_j] : y) {
// write values to file filestream << "\t" << fid_j[i];
auto size = x.size(); }
for (unsigned int i = 0; i < size; i++) { filestream << "\n";
filestream << x[i] << "\t" << y[i] << "\n";
}
} }
}
}
void save_data_to_file(
const std::string &resulttype, const std::string &directory,
const std::vector<double> &x, const std::vector<double> &y,
const std::unordered_map<std::string, double> &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";
}
}
} }

39
src/utils/io.h
View File

@@ -1,29 +1,38 @@
#ifndef RWSIM_IO_H #ifndef RWSIM_IO_H
#define RWSIM_IO_H #define RWSIM_IO_H
#include <unordered_map> #include <filesystem>
#include <map> #include <map>
#include <string> #include <string>
#include <filesystem> #include <unordered_map>
#include <vector> #include <vector>
struct Arguments { struct Arguments {
std::string parameter_file{}; std::string parameter_file{};
bool ste = false; bool ste = false;
bool spectrum = false; bool spectrum = false;
std::string motion_type{}; std::string motion_type{};
std::string distribution_type{}; std::string distribution_type{};
std::unordered_map<std::string, double> optional; std::unordered_map<std::string, double> optional;
}; };
Arguments parse_args(int argc, char* argv[]); Arguments parse_args(int argc, char *argv[]);
std::pair<std::string, double> get_optional_parameter(std::vector<std::string>::const_iterator &it); std::pair<std::string, double>
std::unordered_map<std::string, double> read_parameter(const std::filesystem::path&); get_optional_parameter(std::vector<std::string>::const_iterator &it);
std::unordered_map<std::string, double>
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<std::string, double>&, const std::unordered_map<std::string, double>&); void save_parameter_to_file(const std::string &, const std::string &,
void save_data_to_file(const std::string&, const std::string&, const std::vector<double>&, const std::map<double, std::vector<double>>&, const std::unordered_map<std::string, double>&); const std::unordered_map<std::string, double> &,
void save_data_to_file(const std::string&, const std::string&, const std::vector<double>&, const std::vector<double>&, const std::unordered_map<std::string, double>&); const std::unordered_map<std::string, double> &);
void save_data_to_file(const std::string &, const std::string &,
const std::vector<double> &,
const std::map<double, std::vector<double>> &,
const std::unordered_map<std::string, double> &);
void save_data_to_file(const std::string &, const std::string &,
const std::vector<double> &, const std::vector<double> &,
const std::unordered_map<std::string, double> &);
#endif #endif

87
src/utils/ranges.cpp
View File

@@ -1,57 +1,58 @@
#include <vector>
#include <algorithm> #include <algorithm>
#include <cmath> #include <cmath>
#include <vector>
#include "ranges.h" #include "ranges.h"
std::vector<double> arange(const int size, const double spacing = 1.) {
std::vector<double> out(size);
std::generate(out.begin(), out.end(),
[n = 0, spacing]() mutable { return n++ * spacing; });
std::vector<double> arange(const int size, const double spacing=1.) { return out;
std::vector<double> out(size);
std::generate(out.begin(), out.end(), [n = 0, spacing]() mutable { return n++ * spacing; });
return out;
} }
std::vector<double> linspace(const double start, const double stop, const int steps) { std::vector<double> linspace(const double start, const double stop,
std::vector<double> range; const int steps) {
range.reserve(steps); std::vector<double> range;
range.reserve(steps);
if (steps == 0) { 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<steps; i++) {
range.push_back(start + stepsize * i);
}
return range; return range;
}
if (steps == 1) {
range.push_back(start);
return range;
}
const double stepsize = (stop - start) / (steps - 1);
for (int i = 0; i < steps; i++) {
range.push_back(start + stepsize * i);
}
return range;
} }
std::vector<double> logspace(const double start, const double stop,
const int steps) {
std::vector<double> range;
range.reserve(steps);
std::vector<double> logspace(const double start, const double stop, const int steps) { if (steps == 0) {
std::vector<double> 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<steps; i++) {
range.push_back(std::pow(10, logstart + stepsize * i));
}
return range; 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 < steps; i++) {
range.push_back(std::pow(10, logstart + stepsize * i));
}
return range;
}

62
src/utils/registry.h
View File

@@ -3,48 +3,46 @@
#include <functional> #include <functional>
#include <memory> #include <memory>
#include <stdexcept>
#include <string> #include <string>
#include <unordered_map> #include <unordered_map>
#include <stdexcept>
#include <vector> #include <vector>
template<typename Base> template <typename Base> class Registry {
class Registry {
public: public:
using Creator = std::function<std::unique_ptr<Base>()>; using Creator = std::function<std::unique_ptr<Base>()>;
static std::unordered_map<std::string, Creator>& entries() { static std::unordered_map<std::string, Creator> &entries() {
static std::unordered_map<std::string, Creator> map; static std::unordered_map<std::string, Creator> map;
return map; return map;
} }
static void add(const std::string& name, Creator creator) { static void add(const std::string &name, Creator creator) {
entries()[name] = std::move(creator); entries()[name] = std::move(creator);
} }
static std::unique_ptr<Base> create(const std::string& name) { static std::unique_ptr<Base> create(const std::string &name) {
auto& map = entries(); auto &map = entries();
auto it = map.find(name); auto it = map.find(name);
if (it == map.end()) { if (it == map.end()) {
std::string msg = "Unknown model '" + name + "'. Available: "; std::string msg = "Unknown model '" + name + "'. Available: ";
for (const auto& [key, _] : map) { for (const auto &[key, _] : map) {
msg += key + ", "; msg += key + ", ";
} }
if (!map.empty()) { if (!map.empty()) {
msg.pop_back(); msg.pop_back();
msg.pop_back(); msg.pop_back();
} }
throw std::invalid_argument(msg); throw std::invalid_argument(msg);
}
return it->second();
} }
return it->second();
}
}; };
template<typename Base, typename Derived> template <typename Base, typename Derived> struct AutoRegister {
struct AutoRegister { explicit AutoRegister(const std::string &name) {
explicit AutoRegister(const std::string& name) { Registry<Base>::add(name, []() { return std::make_unique<Derived>(); });
Registry<Base>::add(name, []() { return std::make_unique<Derived>(); }); }
}
}; };
#endif //RWSIM_REGISTRY_H #endif // RWSIM_REGISTRY_H