From d844aac0e88e7f112b58327e9967471f494cfbda Mon Sep 17 00:00:00 2001 From: Dominik Demuth Date: Thu, 28 Nov 2024 14:50:26 +0100 Subject: [PATCH] added flexibility --- .gitignore | 1 + README.md | 2 -- build.sh | 1 + config.txt | 2 +- main.py | 26 ++++++++++++------- python/helpers.py | 20 ++++++++++---- python/ste.py | 49 ++++++++++++++++------------------- src/motions/base.h | 2 ++ src/motions/bimodalangle.cpp | 1 + src/motions/isosmallangle.cpp | 1 + src/motions/random.cpp | 4 ++- src/motions/tetrahedral.cpp | 1 + src/simulation/sims.cpp | 28 ++++++++++++++++---- src/simulation/sims.h | 3 ++- src/utils/io.cpp | 30 +++++++++++++++++++++ src/utils/io.h | 2 ++ 16 files changed, 121 insertions(+), 52 deletions(-) diff --git a/.gitignore b/.gitignore index dfcdc9f..4a97ccf 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,3 @@ /cmake-build-debug/ .idea +/build/ diff --git a/README.md b/README.md index 0b697c2..cfd88b2 100644 --- a/README.md +++ b/README.md @@ -5,8 +5,6 @@ # Build -This - # Running ## Command line diff --git a/build.sh b/build.sh index 292e835..f8916d8 100644 --- a/build.sh +++ b/build.sh @@ -5,3 +5,4 @@ cmake .. cmake --build . cp ./src/rwsim ../rwsim +cd .. diff --git a/config.txt b/config.txt index 6f6891b..c2a1337 100644 --- a/config.txt +++ b/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 diff --git a/main.py b/main.py index 06ad8c3..582ae99 100644 --- a/main.py +++ b/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) - print(conf_file) + conf_file = find_config_file(motion, distribution, variation) - 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) diff --git a/python/helpers.py b/python/helpers.py index eb4145f..42c5786 100644 --- a/python/helpers.py +++ b/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]: diff --git a/python/ste.py b/python/ste.py index 4ba02f5..ddf5a45 100644 --- a/python/ste.py +++ b/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_ss = np.loadtxt(f'sinsin_{varied_string}.dat') + raw_data = np.loadtxt(f'{prefix}_{varied_string}.dat') - t_mix = raw_data_cc[:, 0] - cc_decay = raw_data_cc[:, 1:] - ss_decay = raw_data_ss[:, 1:] + t_mix = raw_data[:, 0] + decay = raw_data[:, 1:] if plot_decays: - fig_cc, ax_cc = plt.subplots() - ax_cc.set_title('Cos-Cos') - ax_cc.semilogx(t_mix, cc_decay, '.') + fig, ax = plt.subplots() + ax.set_title(prefix) + ax.semilogx(t_mix, decay, '.') - fig_ss, ax_ss = plt.subplots() - ax_ss.set_title('Sin-Sin') - ax_ss.semilogx(t_mix, ss_decay, '.') + fig.show() - 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'beta_cc_{varied_string}.dat', beta_cc) - np.savetxt(f'finfty_cc_{varied_string}.dat', finfty_cc) + np.savetxt(f'tau_{prefix}_{varied_string}.dat', tau) + np.savetxt(f'beta_{prefix}_{varied_string}.dat', beta) + np.savetxt(f'finfty_{prefix}_{varied_string}.dat', finfty) - print('Fit Sin-Sin') - 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 + return varied_string, tau, beta, finfty diff --git a/src/motions/base.h b/src/motions/base.h index a4592b9..bd5056f 100644 --- a/src/motions/base.h +++ b/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); diff --git a/src/motions/bimodalangle.cpp b/src/motions/bimodalangle.cpp index 98aaba1..cf0952f 100644 --- a/src/motions/bimodalangle.cpp +++ b/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() { diff --git a/src/motions/isosmallangle.cpp b/src/motions/isosmallangle.cpp index 2475bd6..b633e92 100644 --- a/src/motions/isosmallangle.cpp +++ b/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() { diff --git a/src/motions/random.cpp b/src/motions/random.cpp index f780c08..5139241 100644 --- a/src/motions/random.cpp +++ b/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()); diff --git a/src/motions/tetrahedral.cpp b/src/motions/tetrahedral.cpp index ead6c0c..c9ef5c5 100644 --- a/src/motions/tetrahedral.cpp +++ b/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() { diff --git a/src/simulation/sims.cpp b/src/simulation/sims.cpp index 30ced7a..5994c12 100644 --- a/src/simulation/sims.cpp +++ b/src/simulation/sims.cpp @@ -49,8 +49,9 @@ void run_spectrum( for (int mol_i = 0; mol_i < num_walker; mol_i++){ std::vector traj_time{}; std::vector traj_phase{}; + std::vector 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(num_mix_times); ss_dict[t_evo_i] = std::vector(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 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 traj_time{}; std::vector traj_phase{}; + std::vector 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& out_time, - std::vector& out_phase + std::vector& out_phase, + std::vector& 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); } } diff --git a/src/simulation/sims.h b/src/simulation/sims.h index 1393113..08fbf5a 100644 --- a/src/simulation/sims.h +++ b/src/simulation/sims.h @@ -36,8 +36,9 @@ void run_ste(std::unordered_map& 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& out_time, std::vector& out_phase); +void make_trajectory(Motion& motion, Distribution& dist, double t_max, std::vector& out_time, std::vector& out_phase, std::vector& out_omega); std::chrono::system_clock::time_point printStart(std::unordered_map &optional); void printEnd(std::chrono::system_clock::time_point start); diff --git a/src/utils/io.cpp b/src/utils/io.cpp index b7afbe9..c9875c5 100644 --- a/src/utils/io.cpp +++ b/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& x, + const std::vector& y, + std::unordered_map& 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"; + } + } +} diff --git a/src/utils/io.h b/src/utils/io.h index e26fbf5..525f494 100644 --- a/src/utils/io.h +++ b/src/utils/io.h @@ -21,7 +21,9 @@ std::pair get_optional_parameter(std::vector:: std::unordered_map read_parameter(const std::filesystem::path&); + void save_parameter_to_file(const std::string&, const std::string&, const std::string&, std::unordered_map&, std::unordered_map&); void save_data_to_file(const std::string&, const std::string&, const std::string&, const std::vector&, const std::map>&, std::unordered_map&); +void save_data_to_file(const std::string&, const std::string&, const std::string&, const std::vector&, const std::vector&, std::unordered_map&); #endif