first commit

.gitignore

@@ -0,0 +1,2 @@
+/cmake-build-debug/
+.idea

CMakeLists.txt

@@ -0,0 +1,27 @@
+cmake_minimum_required(VERSION 3.28)
+project(rwsim)
+
+set(CMAKE_CXX_STANDARD 17)
+
+add_executable(rwsim main.cpp
+        functions.cpp
+        functions.h
+        io.cpp
+        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
+        sims.cpp
+        sims.h
+        ranges.cpp
+        ranges.h
+        motions/tetrahedral.cpp
+        motions/tetrahedral.h
+)
+
+target_compile_options(rwsim PUBLIC -Werror -Wall -Wextra -Wconversion -O2)

functions.cpp

@@ -0,0 +1,51 @@
+//
+// Created by dominik on 8/14/24.
+//#
+
+#include <vector>
+#include <chrono>
+#include <iostream>
+
+
+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) {
+    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(
+    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;
+    }
+
+    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;
+    // time_t end_time = std::chrono::system_clock::to_time_t(end);
+    // std::cout << "End   tau = " << tau_i << "s : " << ctime(&end_time);
+    // std::cout << "Duration: " << duration.count() << "s\n" << std::endl;
+
+}

functions.h

@@ -0,0 +1,13 @@
+#ifndef RWSIM_FUNCTIONS_H
+#define RWSIM_FUNCTIONS_H
+
+#include <vector>
+
+
+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

io.cpp

@@ -0,0 +1,88 @@
+//
+// Created by dominik on 8/14/24.
+//
+#include "io.h"
+
+#include <sstream>
+#include <fstream>
+#include <iostream>
+#include <algorithm>
+#include <vector>
+#include <iomanip>
+#include <unordered_map>
+#include <map>
+#include <string>
+
+
+std::unordered_map<std::string, double> parse_arguments(const char *infile) {
+    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;
+
+    // TODO this needs a check for file existence
+
+    while (std::getline(instream, line)) {
+        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);
+    }
+
+    std::cout << "Found parameter\n";
+    for (const auto& [key, value]: parameter) {
+        std::cout << "  " << key << ": " << std::to_string(value) << "\n";
+    }
+    std::cout << std::endl;
+
+    return parameter;
+}
+
+
+void fid_write_out(const std::string& filename, const std::vector<double>& x, const std::vector<double>& y, const double tau, const double t_evo) {
+    auto size = x.size();
+
+    std::ostringstream sfile;
+    sfile << filename << "_";
+    sfile << std::setprecision(6) << std::scientific;
+    sfile << "tau=" << tau << "_tevo=" << t_evo << ".dat";
+    {
+        std::string outfile = sfile.str();
+        std::ofstream fid_file(outfile, std::ios::out);
+        for (unsigned int i = 0; i < size; i++) {
+            fid_file << x[i] << "\t" << y[i] << "\n";
+        }
+    }
+}
+
+void fid_write_out(const std::string& filename, const std::vector<double>& x, const std::map<double, std::vector<double>>& y, const double tau) {
+    auto size = x.size();
+
+    std::ostringstream sfile;
+    sfile << filename << "_";
+    sfile << std::setprecision(6) << std::scientific;
+    sfile << "tau=" << tau << ".dat";
+    {
+        std::string outfile = sfile.str();
+        std::ofstream fid_file(outfile, std::ios::out);
+        fid_file << "#";
+        for (const auto& [t_echo_j, _] : y) {
+            fid_file << t_echo_j << "\t";
+        }
+        fid_file << std::endl;
+
+        for (unsigned int i = 0; i < size; i++) {
+            fid_file << x[i];
+            for (const auto& [_, fid_j] : y) {
+                fid_file << "\t" << fid_j[i];
+            }
+            fid_file << "\n";
+        }
+    }
+}

io.h

@@ -0,0 +1,15 @@
+
+#ifndef RWSIM_IO_H
+#define RWSIM_IO_H
+
+#include <unordered_map>
+#include <map>
+#include <string>
+#include <vector>
+
+std::unordered_map<std::string, double> parse_arguments(const char *);
+
+void fid_write_out(const std::string&, const std::vector<double>&, const std::vector<double>&, double, double);
+void fid_write_out(const std::string&, const std::vector<double>&, const std::map<double, std::vector<double>>&, double tau);
+
+#endif

main.cpp

@@ -0,0 +1,21 @@
+#include <unordered_map>
+#include <random>
+
+#include "io.h"
+#include "sims.h"
+#include "motions/random.h"
+#include "times/delta.h"
+
+
+int main () {
+    std::unordered_map parameter { parse_arguments("config.txt") };
+
+    std::random_device rd;
+    std::mt19937_64 rng(rd());
+
+    auto motion = RandomJump(rng);
+    auto dist = DeltaDistribution(rng);
+
+    run_spectrum(parameter, motion, dist);
+
+}

motions/base.cpp

@@ -0,0 +1,31 @@
+//
+// Created by dominik on 8/12/24.
+//
+
+#include <cmath>
+#include <random>
+
+#include "base.h"
+
+Motion::Motion(const double delta, const double eta, std::mt19937_64& rng) : m_delta(delta), m_eta(eta), m_rng(rng) {
+    m_uni_dist = std::uniform_real_distribution(0., 1.);
+}
+
+Motion::Motion(std::mt19937_64& rng) : 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 * cos(2.*phi));
+}
+
+double 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 omega_q(cos_theta, phi);
+}
+

motions/base.h

@@ -0,0 +1,35 @@
+//
+// Created by dominik on 8/12/24.
+//
+
+#ifndef RWSIM_MOTIONBASE_H
+#define RWSIM_MOTIONBASE_H
+
+
+#include <random>
+
+class Motion {
+public:
+    virtual ~Motion() = default;
+
+    Motion(double, double, std::mt19937_64&);
+    explicit Motion(std::mt19937_64&);
+
+    double draw_position();
+    [[nodiscard]] double omega_q(double, double) const;
+
+    virtual double jump() = 0;
+
+    [[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; }
+
+private:
+    double m_delta{1.};
+    double m_eta{0.};
+    std::mt19937_64& m_rng;
+    std::uniform_real_distribution<double> m_uni_dist;
+};
+
+#endif //RWSIM_MOTIONBASE_H

motions/random.cpp

@@ -0,0 +1,14 @@
+//
+// Created by dominik on 8/12/24.
+//
+
+#include "random.h"
+
+
+RandomJump::RandomJump(const double delta, const double eta, std::mt19937_64 &rng) : Motion(delta, eta, rng) {}
+
+RandomJump::RandomJump(std::mt19937_64 &rng) : Motion(rng) {}
+
+double RandomJump::jump() {
+    return draw_position();
+}

motions/random.h

@@ -0,0 +1,19 @@
+//
+// Created by dominik on 8/12/24.
+//
+
+#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&);
+
+    double jump() override;
+};
+
+#endif //RWSIM_MOTIONRANDOMJUMP_H

motions/tetrahedral.cpp

@@ -0,0 +1,14 @@
+//
+// Created by dominik on 8/16/24.
+//
+#include <random>
+#include "tetrahedral.h"
+
+
+TetrahedralJump::TetrahedralJump(const double delta, const double eta, std::mt19937_64& rng) : Motion(delta, eta, rng) {}
+
+TetrahedralJump::TetrahedralJump(std::mt19937_64& rng) : Motion(rng) {}
+
+double TetrahedralJump::jump() {
+    return draw_position();
+}

motions/tetrahedral.h

@@ -0,0 +1,20 @@
+//
+// Created by dominik on 8/16/24.
+//
+
+#ifndef RWSIM_MOTIONTETRAHEDRAL_H
+#define RWSIM_MOTIONTETRAHEDRAL_H
+
+#include "base.h"
+#include <random>
+
+class TetrahedralJump final : public Motion {
+public:
+    TetrahedralJump(double, double, std::mt19937_64&);
+    explicit TetrahedralJump(std::mt19937_64&);
+
+    double jump() override;
+};
+
+
+#endif //RWSIM_MOTIONTETRAHEDRAL_H

ranges.cpp

@@ -0,0 +1,60 @@
+//
+// Created by dominik on 8/14/24.
+//
+
+
+#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 = log10(start);
+    const double logstop = 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;
+}

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

sims.cpp

@@ -0,0 +1,184 @@
+//
+// Created by dominik on 8/14/24.
+//
+#include <iostream>
+#include <algorithm>
+#include <unordered_map>
+#include <map>
+#include <string>
+#include <vector>
+#include <cmath>
+#include <chrono>
+
+#include "functions.h"
+#include "motions/base.h"
+#include "times/base.h"
+#include "ranges.h"
+#include "sims.h"
+
+#include <bits/fs_fwd.h>
+
+#include "io.h"
+
+
+void run_spectrum(std::unordered_map<std::string, double>& parameter, Motion& motion, Distribution& dist) {
+    const int num_acq = static_cast<int>(parameter["num_acq"]);
+    const int num_walker = static_cast<int>(parameter["num_walker"]);
+
+    const std::vector<double> correlation_times = logspace(parameter["tau_start"], parameter["tau_stop"], static_cast<int>(parameter["tau_steps"]));
+    motion.setDelta(parameter["delta"]);
+    motion.setEta(parameter["eta"]);
+
+    // time axis for all time signals
+    const auto 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"]));
+
+    std::map<double, std::vector<double>> fid_dict;
+    for (auto t_evo_i: echo_times) {
+        fid_dict[t_evo_i] = std::vector<double>(num_acq);
+        std::fill(fid_dict[t_evo_i].begin(), fid_dict[t_evo_i].end(), 0.);
+    }
+
+    const double tmax = *std::max_element(echo_times.begin(), echo_times.end()) * 2 + t_fid.back();
+
+    for (const auto tau_i: correlation_times) {
+        auto start = std::chrono::system_clock::now();
+        auto last_print_out = std::chrono::system_clock::now();
+        time_t start_time = std::chrono::system_clock::to_time_t(start);
+        std::cout << "Start tau = " <<  tau_i << "s : " << ctime(&start_time);
+
+        dist.setTau(tau_i);
+
+        for (auto& [_, fid_j]: fid_dict) {
+            std::fill(fid_j.begin(), fid_j.end(), 0.);
+        }
+
+        // reset array for each correlation time
+        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_tevo = lerp(traj_time, traj_phase, t_echo_j, current_pos);
+
+                // time axis by echo delay to get time in trajectory
+
+                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] += cos(phase_acq - 2 * phase_tevo) / num_walker;
+                }
+                last_print_out = printSteps(last_print_out, start, num_walker, mol_i);
+            }
+        }
+
+        // write fid to files
+        fid_write_out("fid", t_fid, fid_dict, tau_i);
+
+        auto end = std::chrono::system_clock::now();
+
+        std::chrono::duration<float> duration = end - start;
+        time_t end_time = std::chrono::system_clock::to_time_t(end);
+        std::cout << "End   tau = " << tau_i << "s : " << ctime(&end_time);
+        std::cout << "Duration: " << duration.count() << "s\n" << std::endl;
+    }
+}
+
+void run_ste(std::unordered_map<std::string, double>& parameter, Motion& motion, Distribution& dist) {
+    const int num_acq = static_cast<int>(parameter[std::string("num_acq")]);
+    const int num_walker = static_cast<int>(parameter[std::string("num_walker")]);
+
+    const std::vector<double> correlation_times = logspace(parameter["tau_start"], parameter["tau_stop"], static_cast<int>(parameter["tau_steps"]));
+
+    motion.setDelta(parameter["delta"]);
+    motion.setEta(parameter["eta"]);
+
+    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 = linspace(parameter["tevo_start"], parameter["tevo_stop"], static_cast<int>(parameter["tevo_steps"]));
+
+
+    std::map<double, std::vector<double>> fid_dict;
+    for (auto t_evo_i: evolution_times) {
+        fid_dict[t_evo_i] = std::vector<double>(num_acq);
+        std::fill(fid_dict[t_evo_i].begin(), fid_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());
+
+    for (const auto tau_i: correlation_times) {
+        auto start = std::chrono::system_clock::now();
+        time_t start_time = std::chrono::system_clock::to_time_t(start);
+        std::cout << "Start tau = " <<  tau_i << "s : " << ctime(&start_time);
+
+        dist.setTau(tau_i);
+
+        for (auto& [_, fid_j]: fid_dict) {
+            std::fill(fid_j.begin(), fid_j.end(), 0.);
+        }
+
+        // reset array for each correlation time
+        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, fid_j] : fid_dict) {
+                int current_pos = nearest_index(traj_time, t_evo_j, 0);
+                const double phase_tevo = lerp(traj_time, traj_phase, t_evo_j, current_pos);
+
+                // time axis by echo delay to get time in trajectory
+
+                for (int acq_idx = 0; acq_idx < num_acq; acq_idx++) {
+                    const double real_time = mixing_times[acq_idx] + t_evo_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] += cos(phase_acq - 2 * phase_tevo);
+                }
+            }
+        }
+
+        // write fid to files
+        fid_write_out("ste", mixing_times, fid_dict, tau_i);
+
+        auto end = std::chrono::system_clock::now();
+
+        std::chrono::duration<float> duration = end - start;
+        time_t end_time = std::chrono::system_clock::to_time_t(end);
+        std::cout << "End   tau = " << tau_i << "s : " << ctime(&end_time);
+        std::cout << "Duration: " << duration.count() << "s\n" << std::endl;
+    }
+}
+
+
+void make_trajectory(Motion& motion, const 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;
+
+    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;
+
+//        phase += jump(delta, eta, rng) * t;
+
+        out_time.emplace_back(t_passed);
+        out_phase.emplace_back(phase);
+    }
+}

sims.h

@@ -0,0 +1,17 @@
+//
+// Created by dominik on 8/14/24.
+//
+
+#ifndef RWSIM_SIMS_H
+#define RWSIM_SIMS_H
+
+#include <unordered_map>
+#include <string>
+
+#include "motions/base.h"
+#include "times/base.h"
+
+void run_spectrum(std::unordered_map<std::string, double>& parameter, Motion& motion, Distribution& dist);
+void make_trajectory(Motion&, const Distribution&, double, std::vector<double>&, std::vector<double>&);
+
+#endif //RWSIM_SIMS_H

times/base.cpp

@@ -0,0 +1,13 @@
+//
+// Created by dominik on 8/12/24.
+//
+#include "base.h"
+
+Distribution::Distribution(const double tau, std::mt19937_64 &rng) : m_tau(tau), m_rng(rng) {}
+
+Distribution::Distribution(std::mt19937_64 &rng) : m_rng(rng) {}
+
+double Distribution::tau_wait() const {
+    return std::exponential_distribution(1./m_tau)(m_rng);
+}
+

times/base.h

@@ -0,0 +1,28 @@
+//
+// Created by dominik on 8/12/24.
+//
+
+#ifndef RWSIM_TIMESBASE_H
+#define RWSIM_TIMESBASE_H
+
+#include <random>
+
+class Distribution {
+public:
+    virtual ~Distribution() = default;
+
+    Distribution(double, std::mt19937_64&);
+    explicit Distribution(std::mt19937_64&);
+
+    [[nodiscard]] double getTau() const { return m_tau; }
+    void setTau(const double tau) { m_tau = tau;}
+
+    virtual void draw_tau() = 0;
+    [[nodiscard]] double tau_wait() const;
+
+private:
+    double m_tau{1.};
+    std::mt19937_64& m_rng;
+};
+
+#endif //RWSIM_TIMESBASE_H

times/delta.cpp

@@ -0,0 +1,10 @@
+//
+// Created by dominik on 8/12/24.
+//
+#include "delta.h"
+
+DeltaDistribution::DeltaDistribution(const double tau, std::mt19937_64& rng) : Distribution(tau, rng) {}
+DeltaDistribution::DeltaDistribution(std::mt19937_64& rng) : Distribution(rng) {}
+
+void DeltaDistribution::draw_tau() {}
+

times/delta.h

@@ -0,0 +1,18 @@
+//
+// Created by dominik on 8/12/24.
+//
+
+#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 draw_tau() override;
+};
+
+#endif //RWSIM_TIMESDELTA_H