From f1749a0af23b67ee4b71bd005d8c098bd0b426cd Mon Sep 17 00:00:00 2001 From: Dominik Demuth Date: Sun, 10 Nov 2024 15:52:54 +0100 Subject: [PATCH] set motion parameters together via map --- CITATION.cff | 24 +++++++++++++++++ CMakeLists.txt | 2 +- build.sh | 2 +- io.cpp | 7 ++++- io.h | 1 + main.cpp | 18 +++++++------ motions/base.cpp | 6 +++++ motions/base.h | 3 +++ motions/bimodalangle.cpp | 8 ++++++ motions/bimodalangle.h | 1 + motions/isosmallangle.cpp | 9 +++++++ motions/isosmallangle.h | 1 + sims.cpp | 11 +++----- test.py | 55 +++++++++++++++++++++++++++------------ 14 files changed, 114 insertions(+), 34 deletions(-) create mode 100644 CITATION.cff diff --git a/CITATION.cff b/CITATION.cff new file mode 100644 index 0000000..314ad28 --- /dev/null +++ b/CITATION.cff @@ -0,0 +1,24 @@ +# This CITATION.cff file was generated with cffinit. +# Visit https://bit.ly/cffinit to generate yours today! + +cff-version: 1.2.0 +title: NMR Random-Walk simulations +message: >- + If you use this software, please cite it using the + metadata from this file. +type: software +authors: + - given-names: Dominik + family-names: Demuth + email: dominik.demuth@tu-darmstadt.de + affiliation: Technische Universität Darmstadt + orcid: 'https://orcid.org/0000-0003-4648-4875' +repository-code: 'https://gitea.pkm.physik.tu-darmstadt.de/NMRRWSims/cpp' +url: >- + https://gitea.pkm.physik.tu-darmstadt.de/NMRRWSims/cpp/README.md +abstract: >- + C++ classes and functions to simulate solid-state NMR + spectra and stimulated echoes +keywords: + - NMR +license: GPL-3.0 diff --git a/CMakeLists.txt b/CMakeLists.txt index a8a2a72..94498c7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,5 +1,5 @@ cmake_minimum_required(VERSION 3.18) -project(rwsim) +project(RWSim VERSION 1.0) set(CMAKE_CXX_STANDARD 17) diff --git a/build.sh b/build.sh index 934ac2c..e726073 100644 --- a/build.sh +++ b/build.sh @@ -1,5 +1,5 @@ mkdir build -cd build +cd build || exit cmake .. cmake --build . diff --git a/io.cpp b/io.cpp index 3c14828..dbd0156 100644 --- a/io.cpp +++ b/io.cpp @@ -17,7 +17,7 @@ Arguments parse_args(const int argc, char **argv) { - if (argc < 2) { + if (argc < 3) { throw std::runtime_error("Not enough arguments: missing parameter file"); } @@ -70,15 +70,20 @@ std::unordered_map read_parameter(const std::filesystem::pa 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); } + // print result std::cout << "Found parameter\n"; for (const auto& [key, value]: parameter) { std::cout << " " << key << ": " << std::to_string(value) << "\n"; diff --git a/io.h b/io.h index bbb303c..5d36c70 100644 --- a/io.h +++ b/io.h @@ -13,6 +13,7 @@ struct Arguments { bool ste = false; bool spectrum = false; std::string motion_type{}; + std::unordered_map optional; }; Arguments parse_args(int argc, char **argv); diff --git a/main.cpp b/main.cpp index 1ddc634..e1716b8 100644 --- a/main.cpp +++ b/main.cpp @@ -1,24 +1,24 @@ -#include -#include -#include + #include "io.h" #include "sims.h" #include "motions/base.h" -#include "motions/bimodalangle.h" -#include "motions/isosmallangle.h" -#include "motions/random.h" -#include "motions/tetrahedral.h" #include "times/delta.h" #include "times/lognormal.h" +#include +#include +#include + + + int main (const int argc, char *argv[]) { Arguments args; try { args = parse_args(argc, argv); - } catch (std::runtime_error& error) { + } catch (std::invalid_argument& error) { std::cerr << error.what() << std::endl; return 1; } @@ -39,4 +39,6 @@ int main (const int argc, char *argv[]) { if (args.ste) { run_ste(parameter, *motion, dist); } + + return 0; } diff --git a/motions/base.cpp b/motions/base.cpp index 901da3f..957703e 100644 --- a/motions/base.cpp +++ b/motions/base.cpp @@ -60,6 +60,12 @@ Motion* Motion::createFromInput(const std::string& input, std::mt19937_64& rng) throw std::invalid_argument("Invalid input " + input); } +void Motion::setParameters(const std::unordered_map ¶meters) { + m_delta = parameters.at("delta"); + m_eta = parameters.at("eta"); +} + + std::ostream& operator<<(std::ostream& os, const Motion& m) { os << m.getName(); return os; diff --git a/motions/base.h b/motions/base.h index e9bbf5b..d48117e 100644 --- a/motions/base.h +++ b/motions/base.h @@ -7,6 +7,7 @@ #include "coordinates.h" #include +#include class Motion { public: @@ -22,6 +23,8 @@ public: virtual void initialize() = 0; virtual double jump() = 0; + virtual void setParameters(const std::unordered_map&); + [[nodiscard]] double getDelta() const { return m_delta; } void setDelta(const double delta) { m_delta = delta; } [[nodiscard]] double getEta() const { return m_eta; } diff --git a/motions/bimodalangle.cpp b/motions/bimodalangle.cpp index e5ef309..0b95231 100644 --- a/motions/bimodalangle.cpp +++ b/motions/bimodalangle.cpp @@ -23,3 +23,11 @@ double BimodalAngle::jump() { return omega_q(m_prev_pos); } + +void BimodalAngle::setParameters(const std::unordered_map ¶meter) { + Motion::setParameters(parameter); + + m_angle1 = parameter.at("angle1") * M_PI / 180.; + m_angle2 = parameter.at("angle2") * M_PI / 180.; + m_prob = parameter.at("probability1"); +} diff --git a/motions/bimodalangle.h b/motions/bimodalangle.h index 09211f7..fc593e1 100644 --- a/motions/bimodalangle.h +++ b/motions/bimodalangle.h @@ -14,6 +14,7 @@ public: void initialize() override; double jump() override; + void setParameters(const std::unordered_map &) override; protected: double m_angle1{0}; diff --git a/motions/isosmallangle.cpp b/motions/isosmallangle.cpp index 4c2bab4..8eb30ea 100644 --- a/motions/isosmallangle.cpp +++ b/motions/isosmallangle.cpp @@ -3,6 +3,10 @@ // #include "isosmallangle.h" + +#include +#include + #include "coordinates.h" @@ -20,3 +24,8 @@ double SmallAngle::jump() { return omega_q(m_prev_pos); } + +void SmallAngle::setParameters(const std::unordered_map ¶meters) { + m_chi = parameters.at("angle") * M_PI / 180.0; + Motion::setParameters(parameters); +} diff --git a/motions/isosmallangle.h b/motions/isosmallangle.h index 8973669..3018100 100644 --- a/motions/isosmallangle.h +++ b/motions/isosmallangle.h @@ -15,6 +15,7 @@ public: void initialize() override; double jump() override; + void setParameters(const std::unordered_map &) override; private: double m_chi{0}; diff --git a/sims.cpp b/sims.cpp index 892f2bb..7e77849 100644 --- a/sims.cpp +++ b/sims.cpp @@ -37,10 +37,9 @@ void run_spectrum(std::unordered_map& parameter, Motion& mo const double tmax = *std::max_element(echo_times.begin(), echo_times.end()) * 2 + t_fid.back(); // set parameter in distribution and motion model - const double tau = parameter["tau"]; + const double tau = parameter.at("tau"); dist.setTau(tau); - motion.setDelta(parameter["delta"]); - motion.setEta(parameter["eta"]); + motion.setParameters(parameter); const auto start = std::chrono::system_clock::now(); auto last_print_out = std::chrono::system_clock::now(); @@ -92,7 +91,6 @@ void run_ste(std::unordered_map& parameter, Motion& motion, 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; @@ -106,10 +104,9 @@ void run_ste(std::unordered_map& parameter, Motion& motion, 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"]; + const double tau = parameter.at("tau"); dist.setTau(tau); - motion.setDelta(parameter["delta"]); - motion.setEta(parameter["eta"]); + motion.setParameters(parameter); const auto start = std::chrono::system_clock::now(); auto last_print_out = std::chrono::system_clock::now(); diff --git a/test.py b/test.py index f4a3156..5f0cdcb 100644 --- a/test.py +++ b/test.py @@ -15,7 +15,7 @@ def run_sims(taus, ste: bool = True, spectrum: bool = False, exec_file: str = '. if spectrum: arguments += ['--spectrum'] - arguments += ['IsotropicAngle'] + arguments += ['BimodalAngle'] with pathlib.Path(config_file).open('a') as f: f.write(f'tau={tau}\n') @@ -69,7 +69,7 @@ def post_process_spectrum(taus, apod, tpulse): def post_process_ste(taus): - tevo = np.linspace(1e-6, 60e-6, num=121) + tevo = np.linspace(1e-6, 120e-6, num=121) for i, tau in enumerate(taus): try: @@ -83,8 +83,8 @@ def post_process_ste(taus): fig_cc_raw, ax_cc_raw = plt.subplots() fig_ss_raw, ax_ss_raw = plt.subplots() - ax_cc_raw.semilogx(t_mix, raw_data_cc[:, 1:]) - ax_ss_raw.semilogx(t_mix, raw_data_ss[:, 1:]) + ax_cc_raw.semilogx(t_mix, raw_data_cc[:, 1:], '.') + ax_ss_raw.semilogx(t_mix, raw_data_ss[:, 1:], '.') scaled_cc = (raw_data_cc[:, 1:]-raw_data_cc[-1, 1:])/(raw_data_cc[0, 1:]-raw_data_cc[-1, 1:]) scaled_ss = (raw_data_ss[:, 1:]-raw_data_ss[-1, 1:])/(raw_data_ss[0, 1:]-raw_data_ss[-1, 1:]) @@ -101,11 +101,15 @@ def post_process_ste(taus): beta_ss_fit = [] finfty_ss_fit = [] + tau_plus_fit = [] + beta_plus_fit = [] + finfty_plus_fit = [] + for j in range(1, raw_data_cc.shape[1]-1): p0 = [ raw_data_cc[0, 1], t_mix[np.argmin(np.abs(scaled_cc[:, j]-np.exp(-1)))], - 1, + 0.5, 0.1 ] @@ -120,7 +124,7 @@ def post_process_ste(taus): p0 = [ raw_data_cc[0, 1], t_mix[np.argmin(np.abs(scaled_ss[:, j]-np.exp(-1)))], - 1, + 0.5, 0.1 ] @@ -132,24 +136,43 @@ def post_process_ste(taus): beta_ss_fit.append(res[0][2]) finfty_ss_fit.append(res[0][3]) - ax_tau.semilogy(tevo[1:], tau_cc_fit, 'C0-') - ax_tau.axhline(tau) + p0 = [ + raw_data_cc[0, 1], + t_mix[np.argmin(np.abs(scaled_ss[:, j]-np.exp(-1)))], + 1, + 0.1 + ] + + res = curve_fit(ste, t_mix, raw_data_ss[:, j+1] + raw_data_cc[:, j+1], p0, bounds=[(0, 0, 0, 0), (np.inf, np.inf, 1, 1)]) + m0, tauc, beta, finfty = res[0] + # print(f'Sin-Sin-Fit for {tevo[j]}: tau_c = {tauc:.6e}, beta={beta:.4e}, amplitude = {m0: .4e}, f_infty={finfty:.4e}') + + tau_plus_fit.append(res[0][1]) + beta_plus_fit.append(res[0][2]) + finfty_plus_fit.append(res[0][3]) + + # ax_tau.axhline(tau) + + ax_tau.semilogy(tevo[1:], np.array(tau_cc_fit)/tau_cc_fit[0], 'C0-') ax_beta.plot(tevo[1:], beta_cc_fit, 'C0-') ax_finfty.plot(tevo[1:], finfty_cc_fit, 'C0-') - ax_tau.semilogy(tevo[1:], tau_ss_fit, 'C1-') - ax_tau.axhline(tau) + ax_tau.semilogy(tevo[1:], np.array(tau_ss_fit)/tau_ss_fit[0], 'C1-') ax_beta.plot(tevo[1:], beta_ss_fit, 'C1-') ax_finfty.plot(tevo[1:], finfty_ss_fit, 'C1-') + ax_tau.semilogy(tevo[1:], np.array(tau_plus_fit)/tau_plus_fit[0], 'C2-') + ax_beta.plot(tevo[1:], beta_plus_fit, 'C2-') + ax_finfty.plot(tevo[1:], finfty_plus_fit, 'C2-') + plt.show() - np.savetxt('cc_tauc.dat', list(zip(tevo[1:], tau_cc_fit))) - np.savetxt('cc_beta.dat', list(zip(tevo[1:], beta_cc_fit))) - np.savetxt('cc_finfty.dat', list(zip(tevo[1:], finfty_cc_fit))) - np.savetxt('ss_tauc.dat', list(zip(tevo[1:], tau_ss_fit))) - np.savetxt('ss_beta.dat', list(zip(tevo[1:], beta_ss_fit))) - np.savetxt('ss_finfty.dat', list(zip(tevo[1:], finfty_ss_fit))) + # np.savetxt('cc_tauc.dat', list(zip(tevo[1:], tau_cc_fit))) + # np.savetxt('cc_beta.dat', list(zip(tevo[1:], beta_cc_fit))) + # np.savetxt('cc_finfty.dat', list(zip(tevo[1:], finfty_cc_fit))) + # np.savetxt('ss_tauc.dat', list(zip(tevo[1:], tau_ss_fit))) + # np.savetxt('ss_beta.dat', list(zip(tevo[1:], beta_ss_fit))) + # np.savetxt('ss_finfty.dat', list(zip(tevo[1:], finfty_ss_fit)))