added flexibility

.gitignore

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

README.md

@@ -5,8 +5,6 @@
 
 # Build
 
-This 
-
 # Running
 
 ## Command line

build.sh

@@ -5,3 +5,4 @@ cmake ..
 cmake --build .
 
 cp ./src/rwsim ../rwsim
+cd ..

config.txt

@@ -4,7 +4,7 @@ num_walker=20000
 delta=126e3
 eta=0.0
 # Distribution part
-tau=1e-2
+tau=1e-3
 angle1=2
 angle2=30
 probability1=0

main.py

@@ -1,11 +1,14 @@
+import matplotlib.pyplot as plt
+
 from python.ste import *
 from python.helpers import *
 
 # Simulation parameter
 motion = 'IsotropicAngle'
 distribution = 'Delta'
+# parameter = {}
 parameter = {
-    "angle": [3, 10, 30, 109.4],
+    "angle": [10, 109.47],
 }
 
 parameter = prepare_rw_parameter(parameter)
@@ -20,27 +23,30 @@ fig_finfty_cc, ax_finfty_cc = plt.subplots()
 ax_finfty_cc.set_title('f_infty_cc')
 
 fig_tau_ss, ax_tau_ss = plt.subplots()
-ax_tau_ss.set_title('tau_cc')
+ax_tau_ss.set_title('tau_ss')
 
 fig_beta_ss, ax_beta_ss = plt.subplots()
-ax_beta_ss.set_title('beta_cc')
+ax_beta_ss.set_title('beta_ss')
 
 fig_finfty_ss, ax_finfty_ss = plt.subplots()
-ax_finfty_ss.set_title('f_infty_cc')
+ax_finfty_ss.set_title('f_infty_ss')
 
 for variation in parameter:
     print(f"\nRun RW for {motion}/{distribution} with arguments {variation}\n")
-    run_sims(motion, distribution, ste=True, spectrum=True, **variation)
+    run_sims(motion, distribution, ste=True, spectrum=False, **variation)
 
-    conf_file = find_config_file(variation)
+    conf_file = find_config_file(motion, distribution, variation)
-    print(conf_file)
 
-    vary_string, tau_cc, beta_cc, finfty_cc, tau_ss, beta_ss, finfty_ss = fit_and_save_ste(conf_file, plot_decays=False, verbose=False)
+    vary_string, tau_cc, beta_cc, finfty_cc = fit_and_save_ste(conf_file, 'coscos', plot_decays=False, verbose=False)
+    _, tau_ss, beta_ss, finfty_ss = fit_and_save_ste(conf_file, 'sinsin', plot_decays=False, verbose=False)
+    _, tau_2, beta_2, finfty_2 = fit_and_save_ste(conf_file, 'f2', plot_decays=True, verbose=True)
 
-    ax_tau_cc.semilogy(*tau_cc.T, label=vary_string)
+    ax_tau_cc.semilogy(tau_cc[:, 0], tau_cc[:, 1], label=vary_string)
+    ax_tau_cc.axhline(tau_2[:, 1], color='k', linestyle='--')
     ax_beta_cc.plot(*beta_cc.T, label=vary_string)
     ax_finfty_cc.plot(*finfty_cc.T, label=vary_string)
-    ax_tau_ss.semilogy(*tau_ss.T, label=vary_string)
+    ax_tau_ss.semilogy(tau_ss[:, 0], tau_ss[:, 1], label=vary_string)
+    ax_tau_ss.axhline(tau_2[:, 1], color='k', linestyle='--')
     ax_beta_ss.plot(*beta_ss.T, label=vary_string)
     ax_finfty_ss.plot(*finfty_ss.T, label=vary_string)
 

python/helpers.py

@@ -53,13 +53,23 @@ def run_sims(
     subprocess.run(arguments)
 
 
-def find_config_file(var_params: dict) -> pathlib.Path:
+def find_config_file(motion: str, distribution: str, var_params: dict) -> pathlib.Path:
     # TODO handle situation if multiple files fit
-    pattern = re.compile('|'.join(([f'{k}={v:1.6e}' for (k, v) in var_params.items()])).replace('.', '\.').replace('+', '\+'))
+    p_file = None
+    if var_params:
+        var_string = '|'.join(([f'{k}={v:1.6e}' for (k, v) in var_params.items()])).replace('.', '\.').replace('+', '\+')
+        pattern = re.compile(var_string)
+        for p_file in pathlib.Path('.').glob('*_parameter.txt'):
+            if len(re.findall(pattern, str(p_file))) == len(var_params) and re.search(f'{motion}_{distribution}', str(p_file)):
+                return p_file
+        raise ValueError(f'No parameter file found for {motion}, {distribution}, {var_params}')
+    else:
+        for p_file in pathlib.Path('.').glob('*_parameter.txt'):
+            if re.search(f'{motion}_{distribution}', str(p_file)):
+                return p_file
+        raise ValueError(f'No parameter file found for {motion}, {distribution}, {var_params}')
+
 
-    for p_file in pathlib.Path('.').glob('*_parameter.txt'):
-        if len(re.findall(pattern, str(p_file))) == len(var_params):
-            return p_file
 
 
 def read_parameter_file(path: str | pathlib.Path) -> dict[str, float]:

python/ste.py

@@ -24,6 +24,8 @@ def ste_decay(x: np.ndarray, m0: float, t: float, beta: float, finfty: float) ->
 
 def fit_decay(x: np.ndarray, y: np.ndarray, tevo: np.ndarray, verbose: bool = True) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
     num_evo = y.shape[1]
+    if num_evo != tevo.size:
+        tevo = np.arange(num_evo)
 
     tau_fit = np.empty((num_evo, 2))
     tau_fit[:, 0] = tevo
@@ -37,7 +39,7 @@ def fit_decay(x: np.ndarray, y: np.ndarray, tevo: np.ndarray, verbose: bool = Tr
     scaled_y = (y-y[-1, :]) / (y[0, :]-y[-1, :])
 
     for j in range(num_evo):
-        p0 = [scaled_y[0, 1], x[np.argmin(np.abs(scaled_y[:, j]-np.exp(-1)))], 0.5, 0.1]
+        p0 = [y[0, 0], x[np.argmin(np.abs(scaled_y[:, j]-np.exp(-1)))], 0.8, 0.1]
 
         try:
             res = curve_fit(ste_decay, x, y[:, j], p0, bounds=[(0, 0, 0., 0), (np.inf, np.inf, 1, 1)])
@@ -57,7 +59,12 @@ def fit_decay(x: np.ndarray, y: np.ndarray, tevo: np.ndarray, verbose: bool = Tr
     return tau_fit, beta_fit, finfty_fit
 
 
-def fit_and_save_ste(parameter_file: pathlib.Path, plot_decays: bool = True, verbose: bool = True):
+def fit_and_save_ste(
+        parameter_file: pathlib.Path,
+        prefix: str,
+        plot_decays: bool = True,
+        verbose: bool = True,
+) -> tuple[str, np.ndarray, np.ndarray, np.ndarray]:
     # read simulation parameters
     parameter = read_parameter_file(parameter_file)
 
@@ -67,36 +74,24 @@ def fit_and_save_ste(parameter_file: pathlib.Path, plot_decays: bool = True, ver
     # make evolution times
     tevo = np.linspace(parameter['tevo_start'], parameter['tevo_stop'], num=int(parameter['tevo_steps']))
 
-    raw_data_cc = np.loadtxt(f'coscos_{varied_string}.dat')
+    raw_data = np.loadtxt(f'{prefix}_{varied_string}.dat')
-    raw_data_ss = np.loadtxt(f'sinsin_{varied_string}.dat')
 
-    t_mix = raw_data_cc[:, 0]
+    t_mix = raw_data[:, 0]
-    cc_decay = raw_data_cc[:, 1:]
+    decay = raw_data[:, 1:]
-    ss_decay = raw_data_ss[:, 1:]
 
     if plot_decays:
-        fig_cc, ax_cc = plt.subplots()
+        fig, ax = plt.subplots()
-        ax_cc.set_title('Cos-Cos')
+        ax.set_title(prefix)
-        ax_cc.semilogx(t_mix, cc_decay, '.')
+        ax.semilogx(t_mix, decay, '.')
 
-        fig_ss, ax_ss = plt.subplots()
+        fig.show()
-        ax_ss.set_title('Sin-Sin')
-        ax_ss.semilogx(t_mix, ss_decay, '.')
 
-        plt.show()
+    print(f'Fit {prefix}')
+    tau, beta, finfty = fit_decay(t_mix, decay, tevo, verbose=verbose)
 
-    print('Fit Cos-Cos')
-    tau_cc, beta_cc, finfty_cc = fit_decay(t_mix, cc_decay, tevo, verbose=verbose)
 
-    np.savetxt(f'tau_cc_{varied_string}.dat', tau_cc)
+    np.savetxt(f'tau_{prefix}_{varied_string}.dat', tau)
-    np.savetxt(f'beta_cc_{varied_string}.dat', beta_cc)
+    np.savetxt(f'beta_{prefix}_{varied_string}.dat', beta)
-    np.savetxt(f'finfty_cc_{varied_string}.dat', finfty_cc)
+    np.savetxt(f'finfty_{prefix}_{varied_string}.dat', finfty)
 
-    print('Fit Sin-Sin')
+    return varied_string, tau, beta, finfty
-    tau_ss, beta_ss, finfty_ss = fit_decay(t_mix, ss_decay, tevo, verbose=verbose)
-
-    np.savetxt(f'tau_ss_{varied_string}.dat', tau_ss)
-    np.savetxt(f'beta_ss_{varied_string}.dat', beta_ss)
-    np.savetxt(f'finfty_ss_{varied_string}.dat', finfty_ss)
-
-    return varied_string, tau_cc, beta_cc, finfty_cc, tau_ss, beta_ss, finfty_ss

src/motions/base.h

@@ -27,6 +27,7 @@ public:
     [[nodiscard]] double getEta() const { return m_eta; }
     void setEta(const double eta) { m_eta = eta; }
     [[nodiscard]] std::string getName() const { return m_name; }
+    [[nodiscard]] double getInitOmega() const { return m_initial_omega; };
 
     static Motion* createFromInput(const std::string& input, std::mt19937_64& rng);
 
@@ -36,6 +37,7 @@ protected:
     double m_eta{0.};
     std::mt19937_64& m_rng;
     std::uniform_real_distribution<> m_uni_dist;
+    double m_initial_omega{0.};
 };
 
 std::ostream& operator<<(std::ostream& os, const Motion& m);

src/motions/bimodalangle.cpp

@@ -11,6 +11,7 @@ BimodalAngle::BimodalAngle(std::mt19937_64 &rng) : Motion(std::string("BimodalAn
 
 void BimodalAngle::initialize() {
     m_prev_pos = draw_position();
+    m_initial_omega = omega_q(m_prev_pos);
 };
 
 double BimodalAngle::jump() {

src/motions/isosmallangle.cpp

@@ -11,6 +11,7 @@ SmallAngle::SmallAngle(std::mt19937_64 &rng) : Motion(std::string("IsotropicAngl
 
 void SmallAngle::initialize() {
     m_prev_pos = draw_position();
+    m_initial_omega = omega_q(m_prev_pos);
 };
 
 double SmallAngle::jump() {

src/motions/random.cpp

@@ -6,7 +6,9 @@ RandomJump::RandomJump(const double delta, const double eta, std::mt19937_64 &rn
 
 RandomJump::RandomJump(std::mt19937_64 &rng) : Motion(std::string("RandomJump"), rng) {}
 
-void RandomJump::initialize() {}
+void RandomJump::initialize() {
+    m_initial_omega = RandomJump::jump();
+}
 
 double RandomJump::jump() {
     return omega_q(draw_position());

src/motions/tetrahedral.cpp

@@ -19,6 +19,7 @@ void TetrahedralJump::initialize() {
         auto corner_pos = rotate(pos, m_beta, alpha + (i-1) * 2*M_PI/3.);
         m_corners[i] = omega_q(corner_pos);
     }
+    m_initial_omega = TetrahedralJump::jump();
 }
 
 double TetrahedralJump::jump() {

src/simulation/sims.cpp

@@ -49,8 +49,9 @@ void run_spectrum(
     for (int mol_i = 0; mol_i < num_walker; mol_i++){
         std::vector<double> traj_time{};
         std::vector<double> traj_phase{};
+        std::vector<double> traj_omega{};
 
-        make_trajectory(motion, dist, tmax, traj_time, traj_phase);
+        make_trajectory(motion, dist, tmax, traj_time, traj_phase, traj_omega);
 
         for (auto& [t_echo_j, fid_j] : fid_dict) {
             // get phase at echo pulse
@@ -96,8 +97,10 @@ void run_ste(
     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(cc_dict[t_evo_i].begin(), cc_dict[t_evo_i].end(), 0.);
         std::fill(ss_dict[t_evo_i].begin(), ss_dict[t_evo_i].end(), 0.);
     }
+    std::vector<double> f2(num_mix_times);
 
     // 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());
@@ -113,8 +116,16 @@ void run_ste(
     for (int mol_i = 0; mol_i < num_walker; mol_i++){
         std::vector<double> traj_time{};
         std::vector<double> traj_phase{};
+        std::vector<double> traj_omega{};
 
-        make_trajectory(motion, dist, tmax, traj_time, traj_phase);
+        make_trajectory(motion, dist, tmax, traj_time, traj_phase, traj_omega);
+
+        int f2_pos = 0;
+        for (int f2_idx=0; f2_idx < num_mix_times; f2_idx++) {
+            const double t_mix_f2 = mixing_times[f2_idx];
+            f2_pos = nearest_index(traj_time, t_mix_f2, f2_pos);
+            f2[f2_idx] += traj_omega[f2_pos] * motion.getInitOmega() / num_walker;
+        }
 
         for (auto& [t_evo_j, _] : cc_dict) {
             auto& cc_j = cc_dict[t_evo_j];
@@ -149,6 +160,7 @@ void run_ste(
     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);
+    save_data_to_file("f2", motion.getName(), dist.getName(), mixing_times, f2, optional);
 
     printEnd(start);
 }
@@ -159,7 +171,8 @@ void make_trajectory(
     Distribution& dist,
     const double t_max,
     std::vector<double>& out_time,
-    std::vector<double>& out_phase
+    std::vector<double>& out_phase,
+    std::vector<double>& out_omega
     ) {
     // Starting position
     double t_passed = 0;
@@ -168,16 +181,21 @@ void make_trajectory(
     motion.initialize();
     dist.initialize();
 
+    double omega = motion.getInitOmega();
+
     out_time.emplace_back(t_passed);
-    out_phase.emplace_back(0);
+    out_phase.emplace_back(phase);
+    out_omega.emplace_back(omega);
 
     while (t_passed < t_max) {
         const double t = dist.tau_wait();
         t_passed += t;
-        phase += motion.jump() * t;
+        omega = motion.jump();
+        phase += omega * t;
 
         out_time.emplace_back(t_passed);
         out_phase.emplace_back(phase);
+        out_omega.emplace_back(omega);
     }
 }
 

src/simulation/sims.h

@@ -36,8 +36,9 @@ void run_ste(std::unordered_map<std::string, double>& parameter, std::unordered_
  * @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
+ * @param out_omega Vector of omega at jump time
  */
-void make_trajectory(Motion& motion, Distribution& dist, double t_max, std::vector<double>& out_time, std::vector<double>& out_phase);
+void make_trajectory(Motion& motion, Distribution& dist, double t_max, std::vector<double>& out_time, std::vector<double>& out_phase, std::vector<double>& out_omega);
 
 std::chrono::system_clock::time_point printStart(std::unordered_map<std::string, double> &optional);
 void printEnd(std::chrono::system_clock::time_point start);

src/utils/io.cpp

@@ -193,3 +193,33 @@ void save_data_to_file(
         }
     }
 }
+
+void save_data_to_file(
+    const std::string& resulttype,
+    const std::string& motiontype,
+    const std::string& disttype,
+    const std::vector<double>& x,
+    const 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);
+
+        // write values to file
+        auto size = x.size();
+        for (unsigned int i = 0; i < size; i++) {
+            filestream << x[i] << "\t" << y[i] << "\n";
+        }
+    }
+}

src/utils/io.h

@@ -21,7 +21,9 @@ std::pair<std::string, double> get_optional_parameter(std::vector<std::string>::
 
 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>&);
+void save_data_to_file(const std::string&, const std::string&, const std::string&, const std::vector<double>&, const std::vector<double>&, std::unordered_map<std::string, double>&);
 
 #endif