From 1b721978be4ab8a0ca109fbd4476d0d2a456cc3e Mon Sep 17 00:00:00 2001 From: Dominik Demuth Date: Tue, 20 Aug 2024 17:51:49 +0200 Subject: [PATCH] fixed problem in spherical_to_xyz --- functions.cpp | 18 ++++ functions.h | 7 ++ io.cpp | 36 +++++++- io.h | 10 ++- main.cpp | 24 ++++-- motions/base.cpp | 5 +- motions/tetrahedral.cpp | 73 ++++++++-------- motions/tetrahedral.h | 7 ++ sims.cpp | 187 +++++++++++++++++++++------------------- sims.h | 1 + 10 files changed, 228 insertions(+), 140 deletions(-) diff --git a/functions.cpp b/functions.cpp index 13c5fed..5e82827 100644 --- a/functions.cpp +++ b/functions.cpp @@ -3,8 +3,11 @@ //# #include +#include #include #include +#include +#include int nearest_index(const std::vector &x_ref, const double x, int start=0) { @@ -28,6 +31,21 @@ double lerp(const std::vector& x_ref, const std::vector& y_ref, return y_left + dydx * (x - x_left); } +std::array spherical_to_xyz(const double cos_theta, const double phi) { + const double sin_theta = std::sin(std::acos(cos_theta)); + const std::array xyz{sin_theta * std::cos(phi), sin_theta * std::sin(phi), cos_theta}; + + return xyz; +} + +std::pair xyz_to_spherical(const std::array& xyz) { + // std::cout << "length " < printSteps( /* * Prints roughly every 10 seconds how many runs were done and gives a time estimation diff --git a/functions.h b/functions.h index a8329a7..8be3e94 100644 --- a/functions.h +++ b/functions.h @@ -2,12 +2,19 @@ #define RWSIM_FUNCTIONS_H #include +#include +#include +#include int nearest_index(const std::vector&, double, int); double lerp(const std::vector&, const std::vector&, double, int); +std::array spherical_to_xyz(double, double); + +std::pair xyz_to_spherical(const std::array& xyz); + std::chrono::time_point printSteps(std::chrono::time_point, std::chrono::time_point, int, int); #endif diff --git a/io.cpp b/io.cpp index 1dcb038..54c8143 100644 --- a/io.cpp +++ b/io.cpp @@ -15,7 +15,39 @@ #include -std::unordered_map parse_arguments(const std::filesystem::path& infile) { + +Arguments parse_args(const int argc, char **argv) { + if (argc < 2) { + throw std::runtime_error("Not enough arguments: missing parameter file"); + } + + Arguments args; + + for (int i=1; i read_parameter(const std::filesystem::path& infile) { if (!std::filesystem::exists(infile)) { std::cerr << "File " << infile << " does not exist" << std::endl; exit(1); @@ -32,6 +64,8 @@ std::unordered_map parse_arguments(const std::filesystem::p size_t delim_pos; while (std::getline(instream, line)) { + if (line[0] == '#' || line.length() == 1) continue; + line.erase(std::remove(line.begin(), line.end(), ' '), line.end()); delim_pos = line.find('='); key = line.substr(0, delim_pos); diff --git a/io.h b/io.h index ec4f620..50f29ad 100644 --- a/io.h +++ b/io.h @@ -8,7 +8,15 @@ #include #include -std::unordered_map parse_arguments(const std::filesystem::path&); +struct Arguments { + std::string parameter_file{}; + bool ste = false; + bool spectrum = false; +}; + +Arguments parse_args(int argc, char **argv); + +std::unordered_map read_parameter(const std::filesystem::path&); void fid_write_out(const std::string&, const std::vector&, const std::vector&, double, double); void fid_write_out(const std::string&, const std::vector&, const std::map>&, double tau); diff --git a/main.cpp b/main.cpp index 98a94e5..33bc3e3 100644 --- a/main.cpp +++ b/main.cpp @@ -5,25 +5,31 @@ #include "io.h" #include "sims.h" #include "motions/random.h" +#include "motions/tetrahedral.h" #include "times/delta.h" -int main (int argc, char *argv[]) { - - if (argc < 2) { - std::cerr << "Usage: " << argv[0] << " PARAMETER_FILE" << std::endl; +int main (const int argc, char *argv[]) { + Arguments args; + try { + args = parse_args(argc, argv); + } catch (std::runtime_error& error) { + std::cerr << error.what() << std::endl; return 1; } - std::cout << argv[0] << std::endl; - std::unordered_map parameter { parse_arguments(argv[1]) }; + std::unordered_map parameter { read_parameter(args.parameter_file) }; std::random_device rd; std::mt19937_64 rng(rd()); - auto motion = RandomJump(rng); + auto motion = TetrahedralJump(rng); auto dist = DeltaDistribution(rng); - run_spectrum(parameter, motion, dist); - + if (args.spectrum) { + run_spectrum(parameter, motion, dist); + } + if (args.ste) { + run_ste(parameter, motion, dist); + } } diff --git a/motions/base.cpp b/motions/base.cpp index 457cc24..5a86d61 100644 --- a/motions/base.cpp +++ b/motions/base.cpp @@ -9,6 +9,9 @@ #include "base.h" +#include +#include + 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.); @@ -22,7 +25,7 @@ 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)); + return M_PI * m_delta * (3. * cos_theta_square - 1. - m_eta * sin_theta_square * std::cos(2.*phi)); } std::pair Motion::draw_position() { diff --git a/motions/tetrahedral.cpp b/motions/tetrahedral.cpp index c9eb2ee..d6fb36d 100644 --- a/motions/tetrahedral.cpp +++ b/motions/tetrahedral.cpp @@ -4,51 +4,50 @@ #include #include "tetrahedral.h" +#include +#include + +#include "../functions.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) {} void TetrahedralJump::initialize() { - // const auto [cos_theta, phi] = draw_position(); - // - // m_corners[0] = omega_q(cos_theta, phi); - // - // const double alpha = 2. * M_PI * m_uni_dist(m_rng); - // std::cout << alpha << std::endl; + const auto [cos_theta, phi] = draw_position(); + m_corners[0] = omega_q(cos_theta, phi); + const double alpha = 2. * M_PI * m_uni_dist(m_rng); - // v1[0] = sin(phi0)*sin(teta0); - // v2[0] = cos(phi0)*sin(teta0); - // v3[0] = cos(teta0); - // - // abs = v1[0]*v1[0]+v2[0]*v2[0]+v3[0]*v3[0]; // Normalization - // v1[0] = v1[0]/abs; - // v2[0] = v2[0]/abs; - // v3[0] = v3[0]/abs; - // - // /*Calculating the components of the other orientationts. Compare Master Thesis M. Sattig,corrected Version (Share/Abschlussarbeiten)*/ - // norm = sqrt(1 - v3[0]*v3[0])+1E-12; - // normI = 1.0/(norm); - // - // cosgama = cos(delta0); - // singama = sin(delta0); - // v1[1] = v1[0]*cosBETA + sinBETA*normI*(-v1[0]*v3[0]*singama - v2[0]*cosgama); - // v2[1] = v2[0]*cosBETA + sinBETA*normI*(-v2[0]*v3[0]*singama + v1[0]*cosgama); - // v3[1] = v3[0]*cosBETA + sinBETA*norm*singama; - // - // cosgama = cos(gama + delta0); - // singama = sin(gama + delta0); - // v1[2] = v1[0]*cosBETA + sinBETA*normI*(-v1[0]*v3[0]*singama - v2[0]*cosgama); - // v2[2] = v2[0]*cosBETA + sinBETA*normI*(-v2[0]*v3[0]*singama + v1[0]*cosgama); - // v3[2] = v3[0]*cosBETA + sinBETA*norm*singama; - // - // cosgama = cos(2.0*gama + delta0); - // singama = sin(2.0*gama + delta0); - // v1[3] = v1[0]*cosBETA + sinBETA*normI*(-v1[0]*v3[0]*singama - v2[0]*cosgama); - // v2[3] = v2[0]*cosBETA + sinBETA*normI*(-v2[0]*v3[0]*singama + v1[0]*cosgama); - // v3[3] = v3[0]*cosBETA + sinBETA*norm*singama; + const auto vec = spherical_to_xyz(cos_theta, phi); + const double norm = std::sqrt(1 - cos_theta * cos_theta) + 1e-15; + + for (int i = 1; i<4; i++) { + const double cos_alpha = std::cos(alpha + (i-1) * 2*M_PI / 3.); + const double sin_alpha = std::sin(alpha + (i-1) * 2*M_PI / 3.); + std::array rotated_position{}; + if (cos_theta != 1 && cos_theta != -1) { + // std::cout << cos_theta << std::endl; + rotated_position = { + 0*m_cos_beta * vec[0] + m_sin_beta * (-vec[0] * vec[2] * sin_alpha - vec[1] * cos_alpha) / norm, + 0*m_cos_beta * vec[1] + m_sin_beta * (-vec[1] * vec[2] * sin_alpha + vec[0] * cos_alpha) / norm, + m_cos_beta * vec[2] + m_sin_beta * norm * sin_alpha + }; + } else { + rotated_position = { + m_sin_beta * cos_alpha, + m_sin_beta * sin_alpha, + m_cos_beta * cos_theta + }; + } + auto [new_cos_theta, new_phi] = xyz_to_spherical(rotated_position); + m_corners[i] = omega_q(new_cos_theta, new_phi); + } } double TetrahedralJump::jump() { - return 0.; + m_corner_idx += m_chooser(m_rng); + m_corner_idx %= 4; + + return m_corners[m_corner_idx]; } diff --git a/motions/tetrahedral.h b/motions/tetrahedral.h index 78cba62..cae0578 100644 --- a/motions/tetrahedral.h +++ b/motions/tetrahedral.h @@ -20,7 +20,14 @@ public: private: const double m_beta{std::acos(-1/3.)}; + const double m_cos_beta{-1./3.}; + const double m_sin_beta{ 2. * std::sqrt(2.)/3. }; + std::array m_corners{}; + int m_corner_idx{0}; + + std::uniform_int_distribution<> m_chooser{1, 3}; + }; diff --git a/sims.cpp b/sims.cpp index c4196ec..efbffde 100644 --- a/sims.cpp +++ b/sims.cpp @@ -10,143 +10,149 @@ #include #include -#include "functions.h" #include "motions/base.h" #include "times/base.h" +#include "functions.h" #include "ranges.h" #include "sims.h" - -#include - #include "io.h" void run_spectrum(std::unordered_map& parameter, Motion& motion, Distribution& dist) { - const int num_acq = static_cast(parameter["num_acq"]); const int num_walker = static_cast(parameter["num_walker"]); - const std::vector correlation_times = logspace(parameter["tau_start"], parameter["tau_stop"], static_cast(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 int num_acq = static_cast(parameter["num_acq"]); + const std::vector t_fid = arange(num_acq, parameter["dwell_time"]); const std::vector echo_times = linspace(parameter["techo_start"], parameter["techo_stop"], static_cast(parameter["techo_steps"])); + // make timesignal vectors and set them to zero std::map> fid_dict; for (auto t_evo_i: echo_times) { fid_dict[t_evo_i] = std::vector(num_acq); std::fill(fid_dict[t_evo_i].begin(), fid_dict[t_evo_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(); - 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 traj_time{}; - std::vector 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] += std::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 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& parameter, Motion& motion, Distribution& dist) { - const int num_acq = static_cast(parameter[std::string("num_acq")]); - const int num_walker = static_cast(parameter[std::string("num_walker")]); - + // set parameter in distribution and motion model const double tau = parameter["tau"]; dist.setTau(tau); - motion.setDelta(parameter["delta"]); motion.setEta(parameter["eta"]); - const std::vector evolution_times = linspace(parameter["tevo_start"], parameter["tevo_stop"], static_cast(parameter["tevo_steps"])); - const std::vector mixing_times = linspace(parameter["tevo_start"], parameter["tevo_stop"], static_cast(parameter["tevo_steps"])); - - std::map> fid_dict; - for (auto t_evo_i: evolution_times) { - fid_dict[t_evo_i] = std::vector(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()); - const auto start = std::chrono::system_clock::now(); + auto last_print_out = std::chrono::system_clock::now(); const time_t start_time = std::chrono::system_clock::to_time_t(start); - std::cout << "Start tau = " << tau << "s : " << ctime(&start_time); + std::cout << "Start tau = " << tau << "s : " << ctime(&start_time); - for (auto& [_, fid_j]: fid_dict) { - std::fill(fid_j.begin(), fid_j.end(), 0.); - } - - // reset array for each correlation time + // let the walker walk for (int mol_i = 0; mol_i < num_walker; mol_i++){ std::vector traj_time{}; std::vector 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); + 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); - // 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; + 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); + 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 - fid_write_out("ste", mixing_times, fid_dict, tau); + fid_write_out("fid", t_fid, fid_dict, tau); + + const auto end = std::chrono::system_clock::now(); + + std::chrono::duration duration = end - start; + const time_t end_time = std::chrono::system_clock::to_time_t(end); + std::cout << "End tau = " << tau << "s : " << ctime(&end_time); + std::cout << "Duration: " << duration.count() << "s\n" << std::endl; + +} + + +void run_ste(std::unordered_map& parameter, Motion& motion, Distribution& dist) { + const int num_walker = static_cast(parameter[std::string("num_walker")]); + + const int num_mix_times = static_cast(parameter[std::string("tmix_steps")]); + const std::vector evolution_times = linspace(parameter["tevo_start"], parameter["tevo_stop"], static_cast(parameter["tevo_steps"])); + const std::vector mixing_times = logspace(parameter["tmix_start"], parameter["tmix_stop"], num_mix_times); + + + // make ste decay vectors and set them to zero + std::map> cc_dict; + std::map> ss_dict; + for (auto t_evo_i: evolution_times) { + cc_dict[t_evo_i] = std::vector(num_mix_times); + ss_dict[t_evo_i] = std::vector(num_mix_times); + std::fill(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 + const double tau = parameter["tau"]; + dist.setTau(tau); + motion.setDelta(parameter["delta"]); + motion.setEta(parameter["eta"]); + + const auto start = std::chrono::system_clock::now(); + const time_t start_time = std::chrono::system_clock::to_time_t(start); + std::cout << "Start tau = " << tau << "s : " << ctime(&start_time); + + // let the walker walk + for (int mol_i = 0; mol_i < num_walker; mol_i++){ + std::vector traj_time{}; + std::vector 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; + } + } + } + + // write to files + fid_write_out("coscos", mixing_times, cc_dict, tau); + fid_write_out("sinsin", mixing_times, ss_dict, tau); const auto end = std::chrono::system_clock::now(); @@ -171,7 +177,6 @@ void make_trajectory(Motion& motion, const Distribution& dist, const double t_ma while (t_passed < t_max) { const double t = dist.tau_wait(); t_passed += t; - phase += motion.jump() * t; out_time.emplace_back(t_passed); diff --git a/sims.h b/sims.h index 61b8386..ee51500 100644 --- a/sims.h +++ b/sims.h @@ -12,6 +12,7 @@ #include "times/base.h" void run_spectrum(std::unordered_map& parameter, Motion& motion, Distribution& dist); +void run_ste(std::unordered_map& parameter, Motion& motion, Distribution& dist); void make_trajectory(Motion&, const Distribution&, double, std::vector&, std::vector&); #endif //RWSIM_SIMS_H