added flexibility

src/CMakeLists.txt

@@ -0,0 +1,35 @@
+cmake_minimum_required(VERSION 3.18)
+
+set(CMAKE_CXX_STANDARD 17)
+
+
+add_executable(rwsim main.cpp
+        utils/functions.h
+        utils/functions.cpp
+        utils/io.cpp
+        utils/io.h
+        motions/base.cpp
+        motions/base.h
+        motions/random.cpp
+        motions/random.h
+        times/base.cpp
+        times/base.h
+        times/delta.cpp
+        times/delta.h
+        simulation/sims.cpp
+        simulation/sims.h
+        utils/ranges.cpp
+        utils/ranges.h
+        motions/tetrahedral.cpp
+        motions/tetrahedral.h
+        motions/isosmallangle.cpp
+        motions/isosmallangle.h
+        motions/coordinates.cpp
+        motions/coordinates.h
+        motions/bimodalangle.cpp
+        motions/bimodalangle.h
+        times/lognormal.cpp
+        times/lognormal.h
+)
+
+target_compile_options(rwsim PUBLIC -Werror -Wall -Wextra -Wconversion -O2)

src/main.cpp

@@ -0,0 +1,50 @@
+
+#include "utils/io.h"
+#include "simulation/sims.h"
+#include "motions/base.h"
+#include "times/base.h"
+
+
+#include <iostream>
+#include <unordered_map>
+#include <random>
+
+
+
+
+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;
+    }
+
+    // print parameter of simulation to inform user
+    std::cout << "Found parameter\n";
+    for (const auto& [key, value]: parameter) {
+        std::cout << key << ": " << std::to_string(value) << "\n";
+    }
+    std::cout << std::endl;
+
+    std::random_device rd;
+    std::mt19937_64 rng(rd());
+
+    Motion *motion = Motion::createFromInput(args.motion_type, rng);
+    Distribution *dist = Distribution::createFromInput(args.distribution_type, rng);
+    if (args.spectrum) {
+       run_spectrum(parameter, args.optional, *motion, *dist);
+    }
+    if (args.ste) {
+        run_ste(parameter, args.optional, *motion, *dist);
+    }
+
+    return 0;
+}

src/motions/base.cpp

@@ -0,0 +1,64 @@
+
+#include "base.h"
+#include "coordinates.h"
+#include "bimodalangle.h"
+#include "isosmallangle.h"
+#include "random.h"
+#include "tetrahedral.h"
+
+#include <stdexcept>
+
+Motion::Motion(std::string name, const double delta, const double eta, std::mt19937_64& rng) : m_name(std::move(name)), m_delta(delta), m_eta(eta), m_rng(rng) {
+    m_uni_dist = std::uniform_real_distribution(0., 1.);
+}
+
+Motion::Motion(std::string name, std::mt19937_64& rng) : m_name(std::move(name)), m_rng(rng) {
+    m_uni_dist = std::uniform_real_distribution(0., 1.);
+}
+
+double Motion::omega_q(const double cos_theta, const double phi) const {
+    const double cos_theta_square = cos_theta * cos_theta;
+    const double sin_theta_square = 1. - cos_theta_square;
+
+    return M_PI * m_delta * (3. * cos_theta_square - 1. - m_eta * sin_theta_square * std::cos(2.*phi));
+}
+
+double Motion::omega_q(const SphericalPos& pos) const {
+    auto [cos_theta, phi] = pos;
+
+    return omega_q(cos_theta, phi);
+}
+
+SphericalPos Motion::draw_position() {
+    const double cos_theta = 1 - 2 * m_uni_dist(m_rng);
+    const double phi = 2.0 * M_PI * m_uni_dist(m_rng);
+
+    return {cos_theta, phi};
+}
+
+Motion* Motion::createFromInput(const std::string& input, std::mt19937_64& rng) {
+    if (input == "TetrahedralJump")
+        return new TetrahedralJump(rng);
+
+    if (input == "IsotropicAngle")
+        return new SmallAngle(rng);
+
+    if (input == "RandomJump")
+        return new RandomJump(rng);
+
+    if (input == "BimodalAngle")
+        return new BimodalAngle(rng);
+
+    throw std::invalid_argument("Invalid input " + input);
+}
+
+void Motion::setParameters(const std::unordered_map<std::string, double> &parameters) {
+    m_delta = parameters.at("delta");
+    m_eta = parameters.at("eta");
+}
+
+
+std::ostream& operator<<(std::ostream& os, const Motion& m) {
+    os << m.getName();
+    return os;
+}

src/motions/base.h

@@ -0,0 +1,43 @@
+#ifndef RWSIM_MOTIONBASE_H
+#define RWSIM_MOTIONBASE_H
+
+#include "coordinates.h"
+
+#include <random>
+#include <unordered_map>
+
+class Motion {
+public:
+    virtual ~Motion() = default;
+
+    Motion(std::string, double, double, std::mt19937_64&);
+    explicit Motion(std::string, std::mt19937_64&);
+
+    SphericalPos draw_position();
+    [[nodiscard]] double omega_q(double, double) const;
+    [[nodiscard]] double omega_q(const SphericalPos&) const;
+
+    virtual void initialize() = 0;
+    virtual double jump() = 0;
+
+    virtual void setParameters(const std::unordered_map<std::string, double>&);
+
+    [[nodiscard]] double getDelta() const { return m_delta; }
+    void setDelta(const double delta) { m_delta = delta; }
+    [[nodiscard]] double getEta() const { return m_eta; }
+    void setEta(const double eta) { m_eta = eta; }
+    [[nodiscard]] std::string getName() const { return m_name; }
+
+    static Motion* createFromInput(const std::string& input, std::mt19937_64& rng);
+
+protected:
+    std::string m_name{"BaseMotion"};
+    double m_delta{1.};
+    double m_eta{0.};
+    std::mt19937_64& m_rng;
+    std::uniform_real_distribution<> m_uni_dist;
+};
+
+std::ostream& operator<<(std::ostream& os, const Motion& m);
+
+#endif //RWSIM_MOTIONBASE_H

src/motions/bimodalangle.cpp

@@ -0,0 +1,30 @@
+
+#include "bimodalangle.h"
+#include "base.h"
+
+BimodalAngle::BimodalAngle(const double delta, const double eta, const double angle1, const double angle2, const double prob, std::mt19937_64 &rng) :
+    Motion(std::string("BimodalAngle"), delta, eta, rng),
+    m_angle1(angle1 * M_PI / 180.0),
+    m_angle2(angle2 * M_PI / 180.0),
+    m_prob(prob) {};
+BimodalAngle::BimodalAngle(std::mt19937_64 &rng) : Motion(std::string("BimodalAngle"), rng) {}
+
+void BimodalAngle::initialize() {
+    m_prev_pos = draw_position();
+};
+
+double BimodalAngle::jump() {
+    const double angle = m_uni_dist(m_rng) < m_prob ? m_angle1 : m_angle2;
+    const double gamma{2 * M_PI * m_uni_dist(m_rng)};
+    m_prev_pos = rotate(m_prev_pos, angle, gamma);
+
+    return omega_q(m_prev_pos);
+}
+
+void BimodalAngle::setParameters(const std::unordered_map<std::string, double> &parameter) {
+    Motion::setParameters(parameter);
+
+    m_angle1 = parameter.at("angle1") * M_PI / 180.;
+    m_angle2 = parameter.at("angle2") * M_PI / 180.;
+    m_prob = parameter.at("probability1");
+}

src/motions/bimodalangle.h

@@ -0,0 +1,23 @@
+
+#ifndef BIMODALANGLE_H
+#define BIMODALANGLE_H
+
+#include "base.h"
+
+class BimodalAngle : public Motion {
+public:
+    BimodalAngle(double, double, double, double, double, std::mt19937_64& );
+    explicit BimodalAngle(std::mt19937_64&);
+
+    void initialize() override;
+    double jump() override;
+    void setParameters(const std::unordered_map<std::string, double> &) override;
+
+protected:
+    double m_angle1{0};
+    double m_angle2{0};
+    double m_prob{0};
+    SphericalPos m_prev_pos{0., 0.};
+};
+
+#endif //BIMODALANGLE_H

src/motions/coordinates.cpp

@@ -0,0 +1,39 @@
+#include "coordinates.h"
+
+#include <cmath>
+#include <iostream>
+
+SphericalPos rotate(const SphericalPos& old_pos, const double alpha, const double beta) {
+    const double sin_alpha{std::sin(alpha)};
+    const double cos_alpha{std::cos(alpha)};
+
+    const double sin_beta{std::sin(beta)};
+    const double cos_beta{std::cos(beta)};
+
+    const double cos_theta{old_pos.cos_theta};
+
+    if (std::abs(cos_theta) == 1) {
+        return xyz_to_spherical(CartesianPos{cos_alpha * cos_beta, cos_alpha * sin_beta, cos_alpha * cos_theta});
+    }
+
+    const double norm{std::sqrt(1 - cos_theta * cos_theta) + 1e-15};
+
+    auto [x, y , z] = spherical_to_xyz(old_pos);
+
+    const auto new_pos = CartesianPos{
+        cos_alpha * x + sin_alpha * (-x * z * sin_beta - y * cos_beta) / norm,
+        cos_alpha * y + sin_alpha * (-y * z * sin_beta + x * cos_beta) / norm,
+        cos_alpha * z + sin_alpha * norm * sin_beta
+    };
+
+    return xyz_to_spherical(new_pos);
+}
+
+CartesianPos spherical_to_xyz(const SphericalPos& pos) {
+    const double sin_theta = std::sin(std::acos(pos.cos_theta));
+    return {sin_theta * std::cos(pos.phi), sin_theta * std::sin(pos.phi), pos.cos_theta};
+}
+
+SphericalPos xyz_to_spherical(const CartesianPos& pos) {
+    return {pos.z, std::atan2(pos.y, pos.x)};
+}

src/motions/coordinates.h

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

src/motions/isosmallangle.cpp

@@ -0,0 +1,26 @@
+#include "isosmallangle.h"
+#include "coordinates.h"
+
+#include <iostream>
+
+
+
+SmallAngle::SmallAngle(const double delta, const double eta, const double chi, std::mt19937_64 &rng) :
+    Motion(std::string("IsotropicAngle"), delta, eta, rng), m_chi(chi * M_PI / 180.0) {};
+SmallAngle::SmallAngle(std::mt19937_64 &rng) : Motion(std::string("IsotropicAngle"), rng) {}
+
+void SmallAngle::initialize() {
+    m_prev_pos = draw_position();
+};
+
+double SmallAngle::jump() {
+    const double gamma{2 * M_PI * m_uni_dist(m_rng)};
+    m_prev_pos = rotate(m_prev_pos, m_chi, gamma);
+
+    return omega_q(m_prev_pos);
+}
+
+void SmallAngle::setParameters(const std::unordered_map<std::string, double> &parameters) {
+    m_chi = parameters.at("angle") * M_PI / 180.0;
+    Motion::setParameters(parameters);
+}

src/motions/isosmallangle.h

@@ -0,0 +1,22 @@
+
+#ifndef RWSIM_MOTIONISOSMALLANGLE_H
+#define RWSIM_MOTIONISOSMALLANGLE_H
+
+#include "base.h"
+#include "coordinates.h"
+
+class SmallAngle final : public Motion {
+public:
+    SmallAngle(double, double, double, std::mt19937_64& );
+    explicit SmallAngle(std::mt19937_64&);
+
+    void initialize() override;
+    double jump() override;
+    void setParameters(const std::unordered_map<std::string, double> &) override;
+
+private:
+    double m_chi{0};
+    SphericalPos m_prev_pos{0., 0.};
+};
+
+#endif //RWSIM_MOTIONISOSMALLANGLE_H

src/motions/random.cpp

@@ -0,0 +1,13 @@
+
+#include "random.h"
+
+
+RandomJump::RandomJump(const double delta, const double eta, std::mt19937_64 &rng) : Motion(std::string("RandomJump"), delta, eta, rng) {}
+
+RandomJump::RandomJump(std::mt19937_64 &rng) : Motion(std::string("RandomJump"), rng) {}
+
+void RandomJump::initialize() {}
+
+double RandomJump::jump() {
+    return omega_q(draw_position());
+}

src/motions/random.h

@@ -0,0 +1,17 @@
+#ifndef RWSIM_MOTIONRANDOMJUMP_H
+#define RWSIM_MOTIONRANDOMJUMP_H
+
+
+#include "base.h"
+#include <random>
+
+class RandomJump final : public Motion {
+public:
+    RandomJump(double, double, std::mt19937_64&);
+    explicit RandomJump(std::mt19937_64&);
+
+    void initialize() override;
+    double jump() override;
+};
+
+#endif //RWSIM_MOTIONRANDOMJUMP_H

src/motions/tetrahedral.cpp

@@ -0,0 +1,29 @@
+#include "tetrahedral.h"
+
+#include <random>
+
+#include "tetrahedral.h"
+
+TetrahedralJump::TetrahedralJump(const double delta, const double eta, std::mt19937_64& rng) :
+    Motion(std::string{"FourSiteTetrahedral"}, delta, eta, rng) {}
+
+TetrahedralJump::TetrahedralJump(std::mt19937_64& rng) : Motion(std::string{"FourSiteTetrahedral"}, rng) {}
+
+void TetrahedralJump::initialize() {
+    const auto pos = draw_position();
+    m_corners[0] = omega_q(pos);
+
+    const double alpha = 2. * M_PI * m_uni_dist(m_rng);
+
+    for (int i = 1; i<4; i++) {
+        auto corner_pos = rotate(pos, m_beta, alpha + (i-1) * 2*M_PI/3.);
+        m_corners[i] = omega_q(corner_pos);
+    }
+}
+
+double TetrahedralJump::jump() {
+    m_corner_idx += m_chooser(m_rng);
+    m_corner_idx %= 4;
+
+    return m_corners[m_corner_idx];
+}

src/motions/tetrahedral.h

@@ -0,0 +1,28 @@
+#ifndef RWSIM_MOTIONTETRAHEDRAL_H
+#define RWSIM_MOTIONTETRAHEDRAL_H
+
+#include "base.h"
+#include <random>
+#include <cmath>
+#include <array>
+
+class TetrahedralJump final : public Motion {
+public:
+    TetrahedralJump(double, double, std::mt19937_64&);
+    explicit TetrahedralJump(std::mt19937_64&);
+
+    void initialize() override;
+    double jump() override;
+
+private:
+    const double m_beta{std::acos(-1/3.)};
+
+    std::array<double, 4> m_corners{};
+    int m_corner_idx{0};
+
+    std::uniform_int_distribution<> m_chooser{1, 3};
+
+};
+
+
+#endif //RWSIM_MOTIONTETRAHEDRAL_H

src/simulation/sims.cpp

@@ -0,0 +1,207 @@
+#include "sims.h"
+#include "../motions/base.h"
+#include "../times/base.h"
+#include "../utils/functions.h"
+#include "../utils/ranges.h"
+#include "../utils/io.h"
+
+#include <iostream>
+#include <algorithm>
+#include <unordered_map>
+#include <map>
+#include <string>
+#include <vector>
+#include <cmath>
+#include <chrono>
+
+
+void run_spectrum(
+    std::unordered_map<std::string, double>& parameter,
+    std::unordered_map<std::string, double> optional,
+    Motion& motion,
+    Distribution& dist
+    ) {
+    const int num_walker = static_cast<int>(parameter["num_walker"]);
+
+    // time axis for all time signals
+    const int num_acq = static_cast<int>(parameter["num_acq"]);
+    const std::vector<double> t_fid = arange(num_acq, parameter["dwell_time"]);
+    const std::vector<double> echo_times = linspace(parameter["techo_start"], parameter["techo_stop"], static_cast<int>(parameter["techo_steps"]));
+
+    // make timesignal vectors and set them to zero
+    std::map<double, std::vector<double>> fid_dict;
+    for (auto t_echo_i: echo_times) {
+        fid_dict[t_echo_i] = std::vector<double>(num_acq);
+        std::fill(fid_dict[t_echo_i].begin(), fid_dict[t_echo_i].end(), 0.);
+    }
+
+    // calculate min length of a trajectory
+    const double tmax = *std::max_element(echo_times.begin(), echo_times.end()) * 2 + t_fid.back();
+
+    // set parameter in distribution and motion model
+    dist.setParameters(parameter);
+    motion.setParameters(parameter);
+
+    const auto start = printStart(optional);
+    auto last_print_out = std::chrono::system_clock::now();
+
+    // let the walker walk
+    for (int mol_i = 0; mol_i < num_walker; mol_i++){
+        std::vector<double> traj_time{};
+        std::vector<double> traj_phase{};
+
+        make_trajectory(motion, dist, tmax, traj_time, traj_phase);
+
+        for (auto& [t_echo_j, fid_j] : fid_dict) {
+            // get phase at echo pulse
+            int current_pos = nearest_index(traj_time, t_echo_j, 0);
+            const double phase_techo = lerp(traj_time, traj_phase, t_echo_j, current_pos);
+
+
+            for (int acq_idx = 0; acq_idx < num_acq; acq_idx++) {
+                const double real_time = t_fid[acq_idx] + 2 * t_echo_j;
+
+                current_pos = nearest_index(traj_time, real_time, current_pos);
+                const double phase_acq = lerp(traj_time, traj_phase, real_time, current_pos);
+
+                fid_j[acq_idx] += std::cos(phase_acq - 2 * phase_techo) / num_walker;
+            }
+            last_print_out = printSteps(last_print_out, start, num_walker, mol_i);
+        }
+    }
+
+    // write fid to files
+    save_parameter_to_file("timesignal", motion.getName(), dist.getName(), parameter, optional);
+    save_data_to_file("timesignal", motion.getName(), dist.getName(), t_fid, fid_dict, optional);
+
+    printEnd(start);
+}
+
+
+void run_ste(
+    std::unordered_map<std::string, double>& parameter,
+    std::unordered_map<std::string, double> optional,
+    Motion& motion,
+    Distribution& dist
+    ) {
+    const int num_walker = static_cast<int>(parameter[std::string("num_walker")]);
+
+    const int num_mix_times = static_cast<int>(parameter[std::string("tmix_steps")]);
+    const std::vector<double> evolution_times = linspace(parameter["tevo_start"], parameter["tevo_stop"], static_cast<int>(parameter["tevo_steps"]));
+    const std::vector<double> mixing_times = logspace(parameter["tmix_start"], parameter["tmix_stop"], num_mix_times);
+
+    // make ste decay vectors and set them to zero
+    std::map<double, std::vector<double>> cc_dict;
+    std::map<double, std::vector<double>> ss_dict;
+    for (auto t_evo_i: evolution_times) {
+        cc_dict[t_evo_i] = std::vector<double>(num_mix_times);
+        ss_dict[t_evo_i] = std::vector<double>(num_mix_times);
+        std::fill(ss_dict[t_evo_i].begin(), ss_dict[t_evo_i].end(), 0.);
+    }
+
+    // each trajectory must have a duration of at least tmax
+    const double tmax = *std::max_element(evolution_times.begin(), evolution_times.end()) * 2 + *std::max_element(mixing_times.begin(), mixing_times.end());
+
+    // set parameter in distribution and motion model
+    dist.setParameters(parameter);
+    motion.setParameters(parameter);
+
+    const auto start = printStart(optional);
+    auto last_print_out = std::chrono::system_clock::now();
+
+    // let the walker walk
+    for (int mol_i = 0; mol_i < num_walker; mol_i++){
+        std::vector<double> traj_time{};
+        std::vector<double> traj_phase{};
+
+        make_trajectory(motion, dist, tmax, traj_time, traj_phase);
+
+        for (auto& [t_evo_j, _] : cc_dict) {
+            auto& cc_j = cc_dict[t_evo_j];
+            auto& ss_j = ss_dict[t_evo_j];
+
+            // get phase at beginning of mixing time
+            int current_pos = nearest_index(traj_time, t_evo_j, 0);
+            const double dephased = lerp(traj_time, traj_phase, t_evo_j, current_pos);
+            const double cc_tevo = std::cos(dephased);
+            const double ss_tevo = std::sin(dephased);
+
+            for (int mix_idx = 0; mix_idx < num_mix_times; mix_idx++) {
+
+                // get phase at end of mixing time
+                const double time_end_mix = mixing_times[mix_idx] + t_evo_j;
+                current_pos = nearest_index(traj_time, time_end_mix, current_pos);
+                const double phase_mix_end = lerp(traj_time, traj_phase, time_end_mix, current_pos);
+
+                // get phase at echo position
+                const double time_echo = mixing_times[mix_idx] + 2 * t_evo_j;
+                current_pos = nearest_index(traj_time, time_echo, current_pos);
+                const double rephased = lerp(traj_time, traj_phase, time_echo, current_pos) - phase_mix_end;
+
+                cc_j[mix_idx] += cc_tevo * std::cos(rephased) / num_walker;
+                ss_j[mix_idx] += ss_tevo * std::sin(rephased) / num_walker;
+            }
+        }
+        last_print_out = printSteps(last_print_out, start, num_walker, mol_i);
+    }
+
+    // write to files
+    save_parameter_to_file("ste", motion.getName(), dist.getName(), parameter, optional);
+    save_data_to_file("coscos", motion.getName(), dist.getName(), mixing_times, cc_dict, optional);
+    save_data_to_file("sinsin", motion.getName(), dist.getName(), mixing_times, ss_dict, optional);
+
+    printEnd(start);
+}
+
+
+void make_trajectory(
+    Motion& motion,
+    Distribution& dist,
+    const double t_max,
+    std::vector<double>& out_time,
+    std::vector<double>& out_phase
+    ) {
+    // Starting position
+    double t_passed = 0;
+    double phase = 0;
+
+    motion.initialize();
+    dist.initialize();
+
+    out_time.emplace_back(t_passed);
+    out_phase.emplace_back(0);
+
+    while (t_passed < t_max) {
+        const double t = dist.tau_wait();
+        t_passed += t;
+        phase += motion.jump() * t;
+
+        out_time.emplace_back(t_passed);
+        out_phase.emplace_back(phase);
+    }
+}
+
+
+std::chrono::system_clock::time_point printStart(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::cout << "Random walk for ";
+    for (const auto& [key, value] : optional) {
+        std::cout << key << " = " << value << "; ";
+    }
+    std::cout << std::endl;
+    std::cout << "Start: " << ctime(&start_time);
+
+    return start;
+}
+
+
+void printEnd(const std::chrono::system_clock::time_point start) {
+    const auto end = std::chrono::system_clock::now();
+
+    const std::chrono::duration<float> duration = end - start;
+    const time_t end_time = std::chrono::system_clock::to_time_t(end);
+    std::cout << "End:   " << ctime(&end_time);
+    std::cout << "Duration: " << duration.count() << "s\n" << std::endl;
+}

src/simulation/sims.h

@@ -0,0 +1,45 @@
+#ifndef RWSIM_SIMS_H
+#define RWSIM_SIMS_H
+
+#include "../motions/base.h"
+#include "../times/base.h"
+
+#include <unordered_map>
+#include <string>
+#include <chrono>
+
+/**
+ * @brief Run simulation for spectra
+ *
+ * @param parameter Dictionary of parameter for simulation
+ * @param optional Dictionary of parameter set via command line
+ * @param motion Motion model
+ * @param dist Distribution of correlation times
+ */
+void run_spectrum(std::unordered_map<std::string, double>& parameter, std::unordered_map<std::string, double> optional, Motion& motion, Distribution& dist);
+
+/**
+ * @brief Run simulation for stimulated echoes
+ *
+ * @param parameter Dictionary of parameter for simulation
+ * @param optional Dictionary of parameter set via command line
+ * @param motion Motion model
+ * @param dist Distribution of correlation times
+ */
+void run_ste(std::unordered_map<std::string, double>& parameter, std::unordered_map<std::string, double> optional, Motion& motion, Distribution& dist);
+
+/**
+ * @brief Create trajectory of a single walker
+ *
+ * @param motion Motion model
+ * @param dist Distribution of correlation times
+ * @param t_max Double that defines maximum time of trajectory
+ * @param out_time Vector of waiting times
+ * @param out_phase Vector of phase between waiting times
+ */
+void make_trajectory(Motion& motion, Distribution& dist, double t_max, std::vector<double>& out_time, std::vector<double>& out_phase);
+
+std::chrono::system_clock::time_point printStart(std::unordered_map<std::string, double> &optional);
+void printEnd(std::chrono::system_clock::time_point start);
+
+#endif //RWSIM_SIMS_H

src/times/base.cpp

@@ -0,0 +1,27 @@
+#include "base.h"
+#include "delta.h"
+#include "lognormal.h"
+
+#include <stdexcept>
+
+Distribution::Distribution(std::string name, const double tau, std::mt19937_64 &rng) : m_name(std::move(name)), m_tau(tau), m_tau_jump(tau), m_rng(rng) {}
+
+Distribution::Distribution(std::string name, std::mt19937_64 &rng) : m_name(std::move(name)), m_rng(rng) {}
+
+double Distribution::tau_wait() const {
+    return std::exponential_distribution(1./m_tau_jump)(m_rng);
+}
+
+void Distribution::setParameters(const std::unordered_map<std::string, double> &parameters) {
+    m_tau = parameters.at("tau");
+}
+
+Distribution* Distribution::createFromInput(const std::string& input, std::mt19937_64& rng) {
+    if (input == "Delta")
+        return new DeltaDistribution(rng);
+
+    if (input == "LogNormal")
+        return new LogNormalDistribution(rng);
+
+    throw std::invalid_argument("Invalid input " + input);
+}

src/times/base.h

@@ -0,0 +1,33 @@
+#ifndef RWSIM_TIMESBASE_H
+#define RWSIM_TIMESBASE_H
+
+#include <random>
+#include <unordered_map>
+
+class Distribution {
+public:
+    virtual ~Distribution() = default;
+
+    Distribution(std::string, double, std::mt19937_64&);
+    explicit Distribution(std::string, std::mt19937_64&);
+
+    [[nodiscard]] double getTau() const { return m_tau; }
+    void setTau(const double tau) { m_tau = tau; }
+    [[nodiscard]] std::string getName() const { return m_name; };
+
+    virtual void setParameters(const std::unordered_map<std::string, double>&);
+
+    virtual void initialize() = 0;
+    virtual void draw_tau() = 0;
+    [[nodiscard]] double tau_wait() const;
+
+    static Distribution* createFromInput(const std::string& input, std::mt19937_64& rng);
+
+protected:
+    std::string m_name{"BaseDistribution"};
+    double m_tau{1.};
+    double m_tau_jump{1.};
+    std::mt19937_64& m_rng;
+};
+
+#endif //RWSIM_TIMESBASE_H

src/times/delta.cpp

@@ -0,0 +1,11 @@
+#include "delta.h"
+
+DeltaDistribution::DeltaDistribution(const double tau, std::mt19937_64& rng) : Distribution(std::string("Delta"), tau, rng) {}
+DeltaDistribution::DeltaDistribution(std::mt19937_64& rng) : Distribution(std::string("Delta"), rng) {}
+
+void DeltaDistribution::initialize() {
+    m_tau_jump = m_tau;
+}
+
+void DeltaDistribution::draw_tau() {}
+

src/times/delta.h

@@ -0,0 +1,15 @@
+#ifndef RWSIM_TIMESDELTA_H
+#define RWSIM_TIMESDELTA_H
+
+#include "base.h"
+
+class DeltaDistribution final : public Distribution {
+public:
+    DeltaDistribution(double, std::mt19937_64&);
+    explicit DeltaDistribution(std::mt19937_64 &rng);
+
+    void initialize() override;
+    void draw_tau() override;
+};
+
+#endif //RWSIM_TIMESDELTA_H

src/times/lognormal.cpp

@@ -0,0 +1,19 @@
+#include "lognormal.h"
+#include <cmath>
+
+LogNormalDistribution::LogNormalDistribution(const double tau, const double sigma, std::mt19937_64& rng) : Distribution(std::string("LogNormal"), tau, rng), m_sigma(sigma), m_distribution(std::log(tau), sigma) {}
+LogNormalDistribution::LogNormalDistribution(std::mt19937_64& rng) : Distribution(std::string("LogNormal"), rng) {}
+
+void LogNormalDistribution::setParameters(const std::unordered_map<std::string, double> &parameters) {
+    m_sigma = parameters.at("sigma");
+    Distribution::setParameters(parameters);
+}
+
+void LogNormalDistribution::initialize() {
+    m_distribution = std::lognormal_distribution(std::log(m_tau), m_sigma);
+    m_tau_jump = m_distribution(m_rng);
+}
+
+void LogNormalDistribution::draw_tau() {
+    m_tau_jump = m_distribution(m_rng);
+}

src/times/lognormal.h

@@ -0,0 +1,23 @@
+#ifndef LOGNORMAL_H
+#define LOGNORMAL_H
+
+#include "base.h"
+#include <random>
+#include <set>
+
+class LogNormalDistribution final : public Distribution {
+public:
+    LogNormalDistribution(double, double, std::mt19937_64&);
+    explicit LogNormalDistribution(std::mt19937_64 &rng);
+
+    void setParameters(const std::unordered_map<std::string, double> &) override;
+
+    void initialize() override;
+    void draw_tau() override;
+
+private:
+    double m_sigma{1};
+    std::lognormal_distribution<> m_distribution;
+};
+
+#endif //LOGNORMAL_H

src/utils/functions.cpp

@@ -0,0 +1,53 @@
+#include <vector>
+#include <array>
+#include <chrono>
+#include <iostream>
+#include <cmath>
+#include <utility>
+
+
+int nearest_index(const std::vector<double> &x_ref, const double x, int start=0) {
+    while (x > x_ref[start+1]) {
+        start++;
+    }
+    return start;
+}
+
+double lerp(const std::vector<double>& x_ref, const std::vector<double>& y_ref, const double x, const int i) {
+    /*
+     * Linear interpolation between two
+     */
+    const double x_left = x_ref[i];
+    const double y_left = y_ref[i];
+    const double x_right = x_ref[i+1];
+    const double y_right = y_ref[i+1];
+
+    const double dydx = (y_right - y_left) / ( x_right - x_left );
+
+    return y_left + dydx * (x - x_left);
+}
+
+
+std::chrono::time_point<std::chrono::system_clock> printSteps(
+    /*
+     * 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> start,
+    const int total,
+    const int steps
+    ) {
+    const auto now = std::chrono::high_resolution_clock::now();
+
+    if (const std::chrono::duration<float> duration = now - last_print_out; duration.count() < 10.) {
+        return last_print_out;
+    }
+
+    const std::chrono::duration<float> duration = now - start;
+    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";
+
+    return now;
+}

src/utils/functions.h

@@ -0,0 +1,14 @@
+#ifndef RWSIM_FUNCTIONS_H
+#define RWSIM_FUNCTIONS_H
+
+#include <vector>
+#include <chrono>
+
+
+int nearest_index(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>, int, int);
+
+#endif

src/utils/io.cpp

@@ -0,0 +1,195 @@
+#include "io.h"
+
+#include <sstream>
+#include <fstream>
+#include <iostream>
+#include <algorithm>
+#include <complex>
+#include <vector>
+#include <iomanip>
+#include <unordered_map>
+#include <map>
+#include <string>
+#include <filesystem>
+
+
+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);
+
+    return std::make_pair(stripped_arg, stripped_value);
+}
+
+/**
+ * @brief Read and parse arguments coming from the command line
+ * @param argc Number of command-line arguments
+ * @param argv List of command-line arguments
+ * @return Arguments
+ */
+Arguments parse_args(const int argc, char* argv[]) {
+    if (argc < 3) {
+        throw std::runtime_error("Not enough arguments: missing parameter file");
+    }
+
+    Arguments 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;
+        }
+
+        // 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()) {
+        throw std::invalid_argument("Missing argument");
+    }
+
+    return args;
+}
+
+
+std::unordered_map<std::string, double> read_parameter(const std::filesystem::path& infile) {
+    if (!std::filesystem::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::string line;
+    std::string delim = "=";
+    std::string key;
+    std::string value;
+    size_t delim_pos;
+
+    while (std::getline(instream, line)) {
+        // skip comment lines starting with #, and empty lines
+        if (line[0] == '#' || line.length() == 1) continue;
+
+        // strip spaces from line to have always key=value
+        line.erase(std::remove(line.begin(), line.end(), ' '), line.end());
+
+        // split at '=' character and add to map
+        delim_pos = line.find('=');
+        key = line.substr(0, delim_pos);
+        value = line.substr(delim_pos+1);
+        parameter[key] = std::stod(value);
+    }
+
+    return parameter;
+}
+
+
+void save_parameter_to_file(
+    const std::string& resulttype,
+    const std::string& motiontype,
+    const std::string& disttype,
+    std::unordered_map<std::string, double>& parameter,
+    std::unordered_map<std::string, double>& optional
+    ) {
+
+    std::ostringstream parameter_filename;
+    parameter_filename << resulttype << "_" << motiontype << "_" << disttype;
+
+    parameter_filename << std::setprecision(6) << std::scientific;
+    for (const auto& [key, value]: optional) {
+        parameter_filename << "_" << key << "=" << value;
+    }
+    parameter_filename << "_parameter.txt";
+
+    {
+        // write data to file, columns are secondary axis (echo delay, evolution times)
+        std::string datafile = parameter_filename.str();
+        std::ofstream filestream(datafile, std::ios::out);
+
+        for (const auto& [key, value]: parameter) {
+            filestream << key << "=" << value << "\n";
+        }
+    }
+}
+
+
+void save_data_to_file(
+    const std::string& resulttype,
+    const std::string& motiontype,
+    const std::string& disttype,
+    const std::vector<double>& x,
+    const std::map<double, std::vector<double>>& y,
+    std::unordered_map<std::string, double>& optional
+    ) {
+    // make file name
+    std::ostringstream datafile_name;
+    datafile_name << resulttype << "_" << motiontype << "_" << disttype;
+    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";
+        }
+    }
+}

src/utils/io.h

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

src/utils/ranges.cpp

@@ -0,0 +1,55 @@
+#include <vector>
+#include <algorithm>
+#include <cmath>
+
+#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; });
+
+    return out;
+}
+
+std::vector<double> linspace(const double start, const double stop, const int steps) {
+    std::vector<double> range;
+
+    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;
+}
+
+
+std::vector<double> logspace(const double start, const double stop, const int steps) {
+    std::vector<double> range;
+
+    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(pow(10, logstart + stepsize * i));
+    }
+    return range;
+}

src/utils/ranges.h

@@ -0,0 +1,12 @@
+#ifndef RWSIM_RANGES_H
+#define RWSIM_RANGES_H
+
+#include <vector>
+
+std::vector<double> arange(int, double);
+
+std::vector<double> linspace(double, double, int);
+
+std::vector<double> logspace(double, double, int);
+
+#endif