set motion parameters together via map

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

CMakeLists.txt

@@ -1,5 +1,5 @@
 cmake_minimum_required(VERSION 3.18)
-project(rwsim)
+project(RWSim VERSION 1.0)
 
 set(CMAKE_CXX_STANDARD 17)
 

build.sh

@@ -1,5 +1,5 @@
 mkdir build
-cd build
+cd build || exit
 
 cmake ..
 cmake --build .

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<std::string, double> 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";

io.h

@@ -13,6 +13,7 @@ struct Arguments {
     bool ste = false;
     bool spectrum = false;
     std::string motion_type{};
+    std::unordered_map<std::string, double> optional;
 };
 
 Arguments parse_args(int argc, char **argv);

main.cpp

@@ -1,24 +1,24 @@
-#include <iostream>
+
-#include <unordered_map>
-#include <random>
 
 #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 <iostream>
+#include <unordered_map>
+#include <random>
+
+
+
 
 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;
 }

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<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;

motions/base.h

@@ -7,6 +7,7 @@
 
 #include "coordinates.h"
 #include <random>
+#include <unordered_map>
 
 class Motion {
 public:
@@ -22,6 +23,8 @@ public:
     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; }

motions/bimodalangle.cpp

@@ -23,3 +23,11 @@ double BimodalAngle::jump() {
 
     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");
+}

motions/bimodalangle.h

@@ -14,6 +14,7 @@ public:
 
     void initialize() override;
     double jump() override;
+    void setParameters(const std::unordered_map<std::string, double> &) override;
 
 protected:
     double m_angle1{0};

motions/isosmallangle.cpp

@@ -3,6 +3,10 @@
 //
 
 #include "isosmallangle.h"
+
+#include <iostream>
+#include <ostream>
+
 #include "coordinates.h"
 
 
@@ -20,3 +24,8 @@ double SmallAngle::jump() {
 
     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);
+}

motions/isosmallangle.h

@@ -15,6 +15,7 @@ public:
 
     void initialize() override;
     double jump() override;
+    void setParameters(const std::unordered_map<std::string, double> &) override;
 
 private:
     double m_chi{0};

sims.cpp

@@ -37,10 +37,9 @@ void run_spectrum(std::unordered_map<std::string, double>& 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.setParameters(parameter);
-    motion.setEta(parameter["eta"]);
 
     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<std::string, double>& parameter, Motion& motion,
     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;
@@ -106,10 +104,9 @@ void run_ste(std::unordered_map<std::string, double>& 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.setParameters(parameter);
-    motion.setEta(parameter["eta"]);
 
     const auto start = std::chrono::system_clock::now();
     auto last_print_out = std::chrono::system_clock::now();

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_cc_raw.semilogx(t_mix, raw_data_cc[:, 1:], '.')
-        ax_ss_raw.semilogx(t_mix, raw_data_ss[:, 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-')
+            p0 = [
-        ax_tau.axhline(tau)
+                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.semilogy(tevo[1:], np.array(tau_ss_fit)/tau_ss_fit[0], 'C1-')
-        ax_tau.axhline(tau)
         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_tauc.dat', list(zip(tevo[1:], tau_cc_fit)))
-        np.savetxt('cc_beta.dat', list(zip(tevo[1:], beta_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('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_tauc.dat', list(zip(tevo[1:], tau_ss_fit)))
-        np.savetxt('ss_beta.dat', list(zip(tevo[1:], beta_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('ss_finfty.dat', list(zip(tevo[1:], finfty_ss_fit)))