diff --git a/CMakeLists.txt b/CMakeLists.txt index 1352afe..22f090e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,33 +3,6 @@ project(RWSim VERSION 1.0) set(CMAKE_CXX_STANDARD 17) -add_executable(rwsim main.cpp - utils/functions.h - utils/functions.cpp - utils/io.cpp - utils/io.h - motions/base.cpp - motions/base.h - motions/random.cpp - motions/random.h - times/base.cpp - times/base.h - times/delta.cpp - times/delta.h - simulation/sims.cpp - simulation/sims.h - utils/ranges.cpp - utils/ranges.h - motions/tetrahedral.cpp - motions/tetrahedral.h - motions/isosmallangle.cpp - motions/isosmallangle.h - motions/coordinates.cpp - motions/coordinates.h - motions/bimodalangle.cpp - motions/bimodalangle.h - times/lognormal.cpp - times/lognormal.h -) +add_subdirectory(src) target_compile_options(rwsim PUBLIC -Werror -Wall -Wextra -Wconversion -O2) diff --git a/README.md b/README.md index 2cc2c12..0b697c2 100644 --- a/README.md +++ b/README.md @@ -3,24 +3,80 @@ 1. Clone repository with `git clone https://gitea.pkm.physik.tu-darmstadt.de/NMRRWSims/cpp.git` 2. After download, change permissions of build.sh in terminal with `chmod a+x build.sh` -# Setup the simulation +# Build -The workflows leaves yet much to be desired, it will become better (maybe) +This -## Changes in `test.py` +# Running + +## Command line + +To run a random walk simulation Execute `rwsim` in the command line with + +```bash + ./rwsim /path/to/config.txt MotionModel DistributionModel --ste --spectrum -ARG1 1 -ARG2 2 +``` + +It needs three positional arguments: the path to your basic configuration file (see below) and the names of the motional model and of the distribution of correlation times. +Set the optional arguments `--ste` and/or `--spectrum` to create Stimulated Echos or normal echo spectra, respectively. +You can overwrite any parameter given in the configuration file by adding it as optional argument with a numerical value, e.g. `-TAU 1e-3` for a correlation time of 1 ms. + + +## Configuration file `config.txt` + +Change the values of delta, eta, mixing times, echo times,... in this file. If a paramter is defined multiple times, only the last one is used. + +### General parameter + +This list of general parameter are necessary for all simulations: + +| Parameter | Description | +|------------|----------------------------| +| num_walker | Number of trajectories | +| delta | Anisotropy parameter in Hz | +| eta | Asymmetry parameter | + +### Distribution of correlation times + +Two distributions are available: A delta distribution `Delta`, i.e. the same correlation time for every walker, and a log-normal distribution `LogNormal`. + ++ Parameters for delta distribution `Delta` + + | Parameter | Description | + |-----------|-----------------------| + | tau | Jump time in s | + ++ Parameters for log-normal distribution `LogNormal` + + | Parameter | Description | + |-----------|--------------------------------------------| + | tau | Maximum jump time of the distribution in s | + | sigma | Standard deviation of the distribution | + +### Motion models + +Four different jump models are implemented: Isotropic random jump `RandomJump`, isotropic jump with a certain angle, isotropic jump with a bimodal distribution of angles, and a tetrahedral jump `TetrahedralJump`. + ++ Isotropic random jump `RandomJump` does not have additional parameters. ++ Tetrahedral jump `TetrahedralJump` does not have additional parameters. ++ Parameters for isotropic jump by an angle `IsotropicAngle` + + | Parameter | Description | + |-----------|----------------------| + | angle | Jump angle in degree | + ++ Parameters for isotropic jump with bimodal angle distribution `BimodalAngle` + + | Parameter | Description | + |--------------|--------------------------------------| + | angle1 | First jump angle in degree | + | angle2 | Second jump angle in degree | + | probability1 | Probality that jump has angle1 (0-1) | -1. To run spectra, add `-spectrum` in line 14 -2. To run STE, add `-ste` in line 14 -3. Change values for tau, line broadening and pulse length in lines 88-90 (to simulate one tau value, set length of array to 1) -4. Comment / Uncomment `post_process_ste` or `post_process_spectrum`. -## Change of models -1. In `main.cpp`, uncomment the motional model and distribution of correlation times you want to use, comment every other model -## Setting more parameters in `config.txt` -1. Change the values of delta, eta, mixing times, echo times,... in this file. If a paramter is defined multiple times, only the last one is used. # Running diff --git a/build.sh b/build.sh index e726073..292e835 100644 --- a/build.sh +++ b/build.sh @@ -4,7 +4,4 @@ cd build || exit cmake .. cmake --build . -cd .. - -python test.py - +cp ./src/rwsim ../rwsim diff --git a/config.txt b/config.txt index 8b1bb4f..6f6891b 100644 --- a/config.txt +++ b/config.txt @@ -4,24 +4,20 @@ num_walker=20000 delta=126e3 eta=0.0 # Distribution part -# this tau value is overwritten if sim is run with test.py -tau=1e-1 +tau=1e-2 angle1=2 angle2=30 -probability1=0.98 +probability1=0 # Spectrum part -dwell_time=1e-6 +dwell_time=1e-8 num_acq=4096 techo_start=0e-6 techo_stop=40e-6 techo_steps=5 # STE part -tevo_start=1e-6 +tevo_start=2e-6 tevo_stop=120e-6 -tevo_steps=8 +tevo_steps=121 tmix_start=1e-5 tmix_stop=1e1 -tmix_steps=31 -tau=0.01 -tau=0.01 -tau=0.01 +tmix_steps=61 diff --git a/main.py b/main.py new file mode 100644 index 0000000..06ad8c3 --- /dev/null +++ b/main.py @@ -0,0 +1,49 @@ +from python.ste import * +from python.helpers import * + +# Simulation parameter +motion = 'IsotropicAngle' +distribution = 'Delta' +parameter = { + "angle": [3, 10, 30, 109.4], +} + +parameter = prepare_rw_parameter(parameter) + +fig_tau_cc, ax_tau_cc = plt.subplots() +ax_tau_cc.set_title('tau_cc') + +fig_beta_cc, ax_beta_cc = plt.subplots() +ax_beta_cc.set_title('beta_cc') + +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') + +fig_beta_ss, ax_beta_ss = plt.subplots() +ax_beta_ss.set_title('beta_cc') + +fig_finfty_ss, ax_finfty_ss = plt.subplots() +ax_finfty_ss.set_title('f_infty_cc') + +for variation in parameter: + print(f"\nRun RW for {motion}/{distribution} with arguments {variation}\n") + run_sims(motion, distribution, ste=True, spectrum=True, **variation) + + conf_file = find_config_file(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) + + ax_tau_cc.semilogy(*tau_cc.T, label=vary_string) + 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_beta_ss.plot(*beta_ss.T, label=vary_string) + ax_finfty_ss.plot(*finfty_ss.T, label=vary_string) + +for ax in [ax_tau_cc, ax_beta_cc, ax_finfty_cc, ax_tau_ss, ax_beta_ss, ax_finfty_ss]: + ax.legend() +plt.show() \ No newline at end of file diff --git a/python/__init__.py b/python/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/python/helpers.py b/python/helpers.py new file mode 100644 index 0000000..eb4145f --- /dev/null +++ b/python/helpers.py @@ -0,0 +1,77 @@ +from __future__ import annotations + +import pathlib +import re +import subprocess +from itertools import product + +def prepare_rw_parameter(parameter: dict) -> list: + """ + Converts a dictionary of iterables to list of dictionaries + + Example: + If Input is {'a': [1, 2, 3], 'b' = [4, 5]}, output is cartesian product of dictionary values, i.e., + [{'a': 1, 'b': 4}, {'a': 1, 'b': 5}, {'a': 2, 'b': 4}, {'a': 2, 'b': 5}, {'a': 3, 'b': 4}, {'a': 3, 'b': 5}] + + :param parameter: dictionary of list values + :return: list of dictionaries + """ + output = {} + for k, v in parameter.items(): + if isinstance(v, (float, int)): + v = [v] + + output[k] = v + + output = list(dict(zip(parameter.keys(), step)) for step in product(*output.values())) + + return output + + +def run_sims( + motion: str, + distribution: str, + ste: bool = True, + spectrum: bool = False, + exec_file: str = './rwsim', + config_file: str = './config.txt', + **kwargs +) -> None: + + # set positional arguments + arguments = [exec_file, config_file, motion, distribution] + + if ste: + arguments += ['--ste'] + if spectrum: + arguments += ['--spectrum'] + + # add optional parameters that overwrite those given by config file + for k, v in kwargs.items(): + arguments += [f'-{k.upper()}', f'{v}'] + + subprocess.run(arguments) + + +def find_config_file(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('+', '\+')) + + 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]: + path = pathlib.Path(path) + if not path.exists(): + raise ValueError(f"No parameter file found at {path}") + + parameter_dict = {} + with path.open('r') as f: + for line in f.readlines(): + k, v = line.split('=') + parameter_dict[k] = float(v) + + k, v = line.split('=') + return parameter_dict diff --git a/python/spectrum.py b/python/spectrum.py new file mode 100644 index 0000000..6fa64b0 --- /dev/null +++ b/python/spectrum.py @@ -0,0 +1,70 @@ +import numpy +import numpy as np +from matplotlib import pyplot + + + +# parameter for spectrum simulations +lb = 2e3 +pulse_length = 2e-6 + + +def dampening(x: np.ndarray, apod: float) -> np.ndarray: + """ + Calculate additional dampening to account e.g. for field inhomogeneities. + + :param x: Time axis in seconds + :param apod: Dampening factor in 1/seconds + :return: Exponential dampening + """ + + return np.exp(-apod * x) + + +def pulse_attn(freq: np.ndarray, t_pulse: float) -> np.ndarray: + """ + Calculate attenuation of signal to account for finite pulse lengths. + + See Schmitt-Rohr/Spieß, eq. 2.126 for more information. + + :param freq: Frequency axis in Hz + :param t_pulse: Assumed pulse length in s + :return: Attenuation factor. + """ + # cf. Schmitt-Rohr/Spieß eq. 2.126; omega_1 * t_p = pi/2 + pi_half_squared = np.pi**2 / 4 + omega = 2 * np.pi * freq + + numerator = np.sin(np.sqrt(pi_half_squared + omega**2 * t_pulse**2 / 2)) + denominator = np.sqrt(pi_half_squared + omega**2 * t_pulse**2 / 4) + + return np.pi * numerator / denominator / 2 + + +def post_process_spectrum(taus, apod, tpulse): + reduction_factor = np.zeros((taus.size, 5)) # hard-coded t_echo :( + + for i, tau in enumerate(taus): + try: + raw_data = np.loadtxt(f'fid_tau={tau:.6e}.dat') + except OSError: + continue + + t = raw_data[:, 0] + timesignal = raw_data[:, 1:] + + timesignal *= dampening(t, apod)[:, None] + timesignal[0, :] /= 2 + + # FT to spectrum + freq = np.fft.fftshift(np.fft.fftfreq(t.size, d=1e-6)) + spec = np.fft.fftshift(np.fft.fft(timesignal, axis=0), axes=0).real + spec *= pulse_attn(freq, t_pulse=tpulse)[:, None] + + reduction_factor[i, :] = 2*timesignal[0, :] + + plt.plot(freq, spec) + plt.show() + + plt.semilogx(taus, reduction_factor, '.') + plt.show() diff --git a/python/ste.py b/python/ste.py new file mode 100644 index 0000000..4ba02f5 --- /dev/null +++ b/python/ste.py @@ -0,0 +1,102 @@ +import pathlib + +import numpy +import numpy as np +from matplotlib import pyplot as plt +from scipy.optimize import curve_fit + +from python.helpers import read_parameter_file + + +def ste_decay(x: np.ndarray, m0: float, t: float, beta: float, finfty: float) -> np.ndarray: + """ + Calculate stimulated-echo decay. + + :param x: Mixing times in seconds. + :param m0: Amplitude + :param t: Correlation time in seconds + :param beta: Stretching parameter + :param finfty: Final plateau + :return: Stimulated-echo decay + """ + return m0 * ((1-finfty) * np.exp(-(x/t)**beta) + finfty) + + +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] + + tau_fit = np.empty((num_evo, 2)) + tau_fit[:, 0] = tevo + + beta_fit = np.empty((num_evo, 2)) + beta_fit[:, 0] = tevo + + finfty_fit = np.empty((num_evo, 2)) + finfty_fit[:, 0] = tevo + + 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] + + try: + res = curve_fit(ste_decay, x, y[:, j], p0, bounds=[(0, 0, 0., 0), (np.inf, np.inf, 1, 1)]) + except RuntimeError as e: + print(f'Fit {j+1} of {num_evo} failed with {e.args}') + continue + + m0, tauc, beta, finfty = res[0] + + if verbose: + print(f'Fit {j+1} of {num_evo}: tau_c = {tauc:.6e}, beta={beta:.4e}, amplitude = {m0: .4e}, f_infty={finfty:.4e}') + + tau_fit[j, 1] = tauc + beta_fit[j, 1] = beta + finfty_fit[j, 1] = finfty + + return tau_fit, beta_fit, finfty_fit + + +def fit_and_save_ste(parameter_file: pathlib.Path, plot_decays: bool = True, verbose: bool = True): + # read simulation parameters + parameter = read_parameter_file(parameter_file) + + # files have form ste_arg=0.000000e+01_parameter.txt, first remove ste part then parameter.txt to get variables + varied_string = str(parameter_file).partition('_')[-1].rpartition('_')[0] + + # 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') + + t_mix = raw_data_cc[:, 0] + cc_decay = raw_data_cc[:, 1:] + ss_decay = raw_data_ss[:, 1:] + + if plot_decays: + fig_cc, ax_cc = plt.subplots() + ax_cc.set_title('Cos-Cos') + ax_cc.semilogx(t_mix, cc_decay, '.') + + fig_ss, ax_ss = plt.subplots() + ax_ss.set_title('Sin-Sin') + ax_ss.semilogx(t_mix, ss_decay, '.') + + plt.show() + + 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) + + 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 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt new file mode 100644 index 0000000..2736a60 --- /dev/null +++ b/src/CMakeLists.txt @@ -0,0 +1,35 @@ +cmake_minimum_required(VERSION 3.18) + +set(CMAKE_CXX_STANDARD 17) + + +add_executable(rwsim main.cpp + utils/functions.h + utils/functions.cpp + utils/io.cpp + utils/io.h + motions/base.cpp + motions/base.h + motions/random.cpp + motions/random.h + times/base.cpp + times/base.h + times/delta.cpp + times/delta.h + simulation/sims.cpp + simulation/sims.h + utils/ranges.cpp + utils/ranges.h + motions/tetrahedral.cpp + motions/tetrahedral.h + motions/isosmallangle.cpp + motions/isosmallangle.h + motions/coordinates.cpp + motions/coordinates.h + motions/bimodalangle.cpp + motions/bimodalangle.h + times/lognormal.cpp + times/lognormal.h +) + +target_compile_options(rwsim PUBLIC -Werror -Wall -Wextra -Wconversion -O2) diff --git a/main.cpp b/src/main.cpp similarity index 78% rename from main.cpp rename to src/main.cpp index bc7022e..5f782cb 100644 --- a/main.cpp +++ b/src/main.cpp @@ -2,8 +2,7 @@ #include "utils/io.h" #include "simulation/sims.h" #include "motions/base.h" -#include "times/delta.h" -#include "times/lognormal.h" +#include "times/base.h" #include @@ -39,15 +38,12 @@ int main (const int argc, char *argv[]) { std::mt19937_64 rng(rd()); Motion *motion = Motion::createFromInput(args.motion_type, rng); - - auto dist = DeltaDistribution(rng); - // auto dist = LogNormalDistribution(1, 2, rng); // arguments: tau_max, sigma - + Distribution *dist = Distribution::createFromInput(args.distribution_type, rng); if (args.spectrum) { - run_spectrum(parameter, *motion, dist); + run_spectrum(parameter, args.optional, *motion, *dist); } if (args.ste) { - run_ste(parameter, *motion, dist); + run_ste(parameter, args.optional, *motion, *dist); } return 0; diff --git a/motions/base.cpp b/src/motions/base.cpp similarity index 94% rename from motions/base.cpp rename to src/motions/base.cpp index 957703e..4b36f8e 100644 --- a/motions/base.cpp +++ b/src/motions/base.cpp @@ -1,21 +1,13 @@ -// -// Created by dominik on 8/12/24. -// - -#include -#include #include "base.h" - -#include -#include - #include "coordinates.h" #include "bimodalangle.h" #include "isosmallangle.h" #include "random.h" #include "tetrahedral.h" +#include + Motion::Motion(std::string name, const double delta, const double eta, std::mt19937_64& rng) : m_name(std::move(name)), m_delta(delta), m_eta(eta), m_rng(rng) { m_uni_dist = std::uniform_real_distribution(0., 1.); } @@ -69,4 +61,4 @@ void Motion::setParameters(const std::unordered_map ¶me std::ostream& operator<<(std::ostream& os, const Motion& m) { os << m.getName(); return os; -} +} \ No newline at end of file diff --git a/motions/base.h b/src/motions/base.h similarity index 93% rename from motions/base.h rename to src/motions/base.h index d48117e..a4592b9 100644 --- a/motions/base.h +++ b/src/motions/base.h @@ -1,11 +1,8 @@ -// -// Created by dominik on 8/12/24. -// - #ifndef RWSIM_MOTIONBASE_H #define RWSIM_MOTIONBASE_H #include "coordinates.h" + #include #include @@ -34,7 +31,7 @@ public: static Motion* createFromInput(const std::string& input, std::mt19937_64& rng); protected: - std::string m_name{"Base motion"}; + std::string m_name{"BaseMotion"}; double m_delta{1.}; double m_eta{0.}; std::mt19937_64& m_rng; diff --git a/motions/bimodalangle.cpp b/src/motions/bimodalangle.cpp similarity index 87% rename from motions/bimodalangle.cpp rename to src/motions/bimodalangle.cpp index 0b95231..98aaba1 100644 --- a/motions/bimodalangle.cpp +++ b/src/motions/bimodalangle.cpp @@ -1,16 +1,13 @@ -// -// Created by dominik on 8/23/24. -// #include "bimodalangle.h" #include "base.h" BimodalAngle::BimodalAngle(const double delta, const double eta, const double angle1, const double angle2, const double prob, std::mt19937_64 &rng) : - Motion(std::string("Bimodal Angle Jump"), delta, eta, rng), + Motion(std::string("BimodalAngle"), delta, eta, rng), m_angle1(angle1 * M_PI / 180.0), m_angle2(angle2 * M_PI / 180.0), m_prob(prob) {}; -BimodalAngle::BimodalAngle(std::mt19937_64 &rng) : Motion(std::string("Bimodal Angle Jump"), rng) {} +BimodalAngle::BimodalAngle(std::mt19937_64 &rng) : Motion(std::string("BimodalAngle"), rng) {} void BimodalAngle::initialize() { m_prev_pos = draw_position(); diff --git a/motions/bimodalangle.h b/src/motions/bimodalangle.h similarity index 92% rename from motions/bimodalangle.h rename to src/motions/bimodalangle.h index fc593e1..a46afb7 100644 --- a/motions/bimodalangle.h +++ b/src/motions/bimodalangle.h @@ -1,6 +1,3 @@ -// -// Created by dominik on 8/23/24. -// #ifndef BIMODALANGLE_H #define BIMODALANGLE_H diff --git a/motions/coordinates.cpp b/src/motions/coordinates.cpp similarity index 95% rename from motions/coordinates.cpp rename to src/motions/coordinates.cpp index db6629d..f7dbdf4 100644 --- a/motions/coordinates.cpp +++ b/src/motions/coordinates.cpp @@ -1,11 +1,7 @@ -// -// Created by dominik on 8/22/24. -// - #include "coordinates.h" + #include #include -#include SphericalPos rotate(const SphericalPos& old_pos, const double alpha, const double beta) { const double sin_alpha{std::sin(alpha)}; diff --git a/motions/coordinates.h b/src/motions/coordinates.h similarity index 90% rename from motions/coordinates.h rename to src/motions/coordinates.h index 18d9649..56b64de 100644 --- a/motions/coordinates.h +++ b/src/motions/coordinates.h @@ -1,6 +1,3 @@ -// -// Created by dominik on 8/22/24. -// #ifndef COORDINATES_H #define COORDINATES_H diff --git a/motions/isosmallangle.cpp b/src/motions/isosmallangle.cpp similarity index 77% rename from motions/isosmallangle.cpp rename to src/motions/isosmallangle.cpp index 8eb30ea..2475bd6 100644 --- a/motions/isosmallangle.cpp +++ b/src/motions/isosmallangle.cpp @@ -1,18 +1,13 @@ -// -// Created by dominik on 8/21/24. -// - #include "isosmallangle.h" +#include "coordinates.h" #include -#include -#include "coordinates.h" SmallAngle::SmallAngle(const double delta, const double eta, const double chi, std::mt19937_64 &rng) : - Motion(std::string("Isotropic Angle Jump"), delta, eta, rng), m_chi(chi * M_PI / 180.0) {}; -SmallAngle::SmallAngle(std::mt19937_64 &rng) : Motion(std::string("Isotropic Angle Jump"), rng) {} + Motion(std::string("IsotropicAngle"), delta, eta, rng), m_chi(chi * M_PI / 180.0) {}; +SmallAngle::SmallAngle(std::mt19937_64 &rng) : Motion(std::string("IsotropicAngle"), rng) {} void SmallAngle::initialize() { m_prev_pos = draw_position(); diff --git a/motions/isosmallangle.h b/src/motions/isosmallangle.h similarity index 92% rename from motions/isosmallangle.h rename to src/motions/isosmallangle.h index 3018100..0f26130 100644 --- a/motions/isosmallangle.h +++ b/src/motions/isosmallangle.h @@ -1,6 +1,3 @@ -// -// Created by dominik on 8/21/24. -// #ifndef RWSIM_MOTIONISOSMALLANGLE_H #define RWSIM_MOTIONISOSMALLANGLE_H diff --git a/motions/random.cpp b/src/motions/random.cpp similarity index 64% rename from motions/random.cpp rename to src/motions/random.cpp index 6457874..f780c08 100644 --- a/motions/random.cpp +++ b/src/motions/random.cpp @@ -1,13 +1,10 @@ -// -// Created by dominik on 8/12/24. -// #include "random.h" -RandomJump::RandomJump(const double delta, const double eta, std::mt19937_64 &rng) : Motion(std::string("Random Jump"), delta, eta, rng) {} +RandomJump::RandomJump(const double delta, const double eta, std::mt19937_64 &rng) : Motion(std::string("RandomJump"), delta, eta, rng) {} -RandomJump::RandomJump(std::mt19937_64 &rng) : Motion(std::string("Random Jump"), rng) {} +RandomJump::RandomJump(std::mt19937_64 &rng) : Motion(std::string("RandomJump"), rng) {} void RandomJump::initialize() {} diff --git a/motions/random.h b/src/motions/random.h similarity index 89% rename from motions/random.h rename to src/motions/random.h index 17e003e..eca43e0 100644 --- a/motions/random.h +++ b/src/motions/random.h @@ -1,7 +1,3 @@ -// -// Created by dominik on 8/12/24. -// - #ifndef RWSIM_MOTIONRANDOMJUMP_H #define RWSIM_MOTIONRANDOMJUMP_H diff --git a/motions/tetrahedral.cpp b/src/motions/tetrahedral.cpp similarity index 82% rename from motions/tetrahedral.cpp rename to src/motions/tetrahedral.cpp index bd21d37..ead6c0c 100644 --- a/motions/tetrahedral.cpp +++ b/src/motions/tetrahedral.cpp @@ -1,14 +1,13 @@ -// -// Created by dominik on 8/16/24. -// +#include "tetrahedral.h" + #include #include "tetrahedral.h" TetrahedralJump::TetrahedralJump(const double delta, const double eta, std::mt19937_64& rng) : - Motion(std::string{"TetrahedralJump"}, delta, eta, rng) {} + Motion(std::string{"FourSiteTetrahedral"}, delta, eta, rng) {} -TetrahedralJump::TetrahedralJump(std::mt19937_64& rng) : Motion(std::string{"TetrahedralJump"}, rng) {} +TetrahedralJump::TetrahedralJump(std::mt19937_64& rng) : Motion(std::string{"FourSiteTetrahedral"}, rng) {} void TetrahedralJump::initialize() { const auto pos = draw_position(); diff --git a/motions/tetrahedral.h b/src/motions/tetrahedral.h similarity index 93% rename from motions/tetrahedral.h rename to src/motions/tetrahedral.h index ec564cb..1872b0e 100644 --- a/motions/tetrahedral.h +++ b/src/motions/tetrahedral.h @@ -1,7 +1,3 @@ -// -// Created by dominik on 8/16/24. -// - #ifndef RWSIM_MOTIONTETRAHEDRAL_H #define RWSIM_MOTIONTETRAHEDRAL_H diff --git a/simulation/sims.cpp b/src/simulation/sims.cpp similarity index 76% rename from simulation/sims.cpp rename to src/simulation/sims.cpp index cc2b190..30ced7a 100644 --- a/simulation/sims.cpp +++ b/src/simulation/sims.cpp @@ -15,8 +15,12 @@ #include - -void run_spectrum(std::unordered_map& parameter, Motion& motion, Distribution& dist) { +void run_spectrum( + std::unordered_map& parameter, + std::unordered_map optional, + Motion& motion, + Distribution& dist + ) { const int num_walker = static_cast(parameter["num_walker"]); // time axis for all time signals @@ -26,23 +30,20 @@ void run_spectrum(std::unordered_map& parameter, Motion& mo // 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.); + for (auto t_echo_i: echo_times) { + fid_dict[t_echo_i] = std::vector(num_acq); + std::fill(fid_dict[t_echo_i].begin(), fid_dict[t_echo_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(); // set parameter in distribution and motion model - const double tau = parameter.at("tau"); - dist.setTau(tau); + dist.setParameters(parameter); motion.setParameters(parameter); - const auto start = std::chrono::system_clock::now(); + const auto start = printStart(optional); 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); // let the walker walk for (int mol_i = 0; mol_i < num_walker; mol_i++){ @@ -70,19 +71,19 @@ void run_spectrum(std::unordered_map& parameter, Motion& mo } // write fid to files - 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; + save_parameter_to_file("timesignal", motion.getName(), dist.getName(), parameter, optional); + save_data_to_file("timesignal", motion.getName(), dist.getName(), t_fid, fid_dict, optional); + printEnd(start); } -void run_ste(std::unordered_map& parameter, Motion& motion, Distribution& dist) { +void run_ste( + std::unordered_map& parameter, + std::unordered_map optional, + 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")]); @@ -102,14 +103,11 @@ 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.at("tau"); - dist.setTau(tau); + dist.setParameters(parameter); motion.setParameters(parameter); - const auto start = std::chrono::system_clock::now(); + const auto start = printStart(optional); 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); // let the walker walk for (int mol_i = 0; mol_i < num_walker; mol_i++){ @@ -148,20 +146,21 @@ void run_ste(std::unordered_map& parameter, Motion& motion, } // 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(); - - const 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; + 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); + printEnd(start); } -void make_trajectory(Motion& motion, Distribution& dist, const double t_max, std::vector& out_time, std::vector& out_phase) { +void make_trajectory( + Motion& motion, + Distribution& dist, + const double t_max, + std::vector& out_time, + std::vector& out_phase + ) { // Starting position double t_passed = 0; double phase = 0; @@ -181,3 +180,28 @@ void make_trajectory(Motion& motion, Distribution& dist, const double t_max, std out_phase.emplace_back(phase); } } + + +std::chrono::system_clock::time_point printStart(std::unordered_map &optional) { + const auto start = std::chrono::system_clock::now(); + const time_t start_time = std::chrono::system_clock::to_time_t(start); + + std::cout << "Random walk for "; + for (const auto& [key, value] : optional) { + std::cout << key << " = " << value << "; "; + } + std::cout << std::endl; + std::cout << "Start: " << ctime(&start_time); + + return start; +} + + +void printEnd(const std::chrono::system_clock::time_point start) { + const auto end = std::chrono::system_clock::now(); + + const std::chrono::duration duration = end - start; + const time_t end_time = std::chrono::system_clock::to_time_t(end); + std::cout << "End: " << ctime(&end_time); + std::cout << "Duration: " << duration.count() << "s\n" << std::endl; +} diff --git a/simulation/sims.h b/src/simulation/sims.h similarity index 69% rename from simulation/sims.h rename to src/simulation/sims.h index e98e047..1393113 100644 --- a/simulation/sims.h +++ b/src/simulation/sims.h @@ -1,7 +1,3 @@ -// -// Created by dominik on 8/14/24. -// - #ifndef RWSIM_SIMS_H #define RWSIM_SIMS_H @@ -10,24 +6,27 @@ #include #include +#include /** * @brief Run simulation for spectra * * @param parameter Dictionary of parameter for simulation + * @param optional Dictionary of parameter set via command line * @param motion Motion model * @param dist Distribution of correlation times */ -void run_spectrum(std::unordered_map& parameter, Motion& motion, Distribution& dist); +void run_spectrum(std::unordered_map& parameter, std::unordered_map optional, Motion& motion, Distribution& dist); /** * @brief Run simulation for stimulated echoes * * @param parameter Dictionary of parameter for simulation + * @param optional Dictionary of parameter set via command line * @param motion Motion model * @param dist Distribution of correlation times */ -void run_ste(std::unordered_map& parameter, Motion& motion, Distribution& dist); +void run_ste(std::unordered_map& parameter, std::unordered_map optional, Motion& motion, Distribution& dist); /** * @brief Create trajectory of a single walker @@ -40,4 +39,7 @@ void run_ste(std::unordered_map& parameter, Motion& motion, */ void make_trajectory(Motion& motion, Distribution& dist, double t_max, std::vector& out_time, std::vector& out_phase); +std::chrono::system_clock::time_point printStart(std::unordered_map &optional); +void printEnd(std::chrono::system_clock::time_point start); + #endif //RWSIM_SIMS_H diff --git a/src/times/base.cpp b/src/times/base.cpp new file mode 100644 index 0000000..bd6f993 --- /dev/null +++ b/src/times/base.cpp @@ -0,0 +1,27 @@ +#include "base.h" +#include "delta.h" +#include "lognormal.h" + +#include + +Distribution::Distribution(std::string name, const double tau, std::mt19937_64 &rng) : m_name(std::move(name)), m_tau(tau), m_tau_jump(tau), m_rng(rng) {} + +Distribution::Distribution(std::string name, std::mt19937_64 &rng) : m_name(std::move(name)), m_rng(rng) {} + +double Distribution::tau_wait() const { + return std::exponential_distribution(1./m_tau_jump)(m_rng); +} + +void Distribution::setParameters(const std::unordered_map ¶meters) { + m_tau = parameters.at("tau"); +} + +Distribution* Distribution::createFromInput(const std::string& input, std::mt19937_64& rng) { + if (input == "Delta") + return new DeltaDistribution(rng); + + if (input == "LogNormal") + return new LogNormalDistribution(rng); + + throw std::invalid_argument("Invalid input " + input); +} \ No newline at end of file diff --git a/src/times/base.h b/src/times/base.h new file mode 100644 index 0000000..4fa08bb --- /dev/null +++ b/src/times/base.h @@ -0,0 +1,33 @@ +#ifndef RWSIM_TIMESBASE_H +#define RWSIM_TIMESBASE_H + +#include +#include + +class Distribution { +public: + virtual ~Distribution() = default; + + Distribution(std::string, double, std::mt19937_64&); + explicit Distribution(std::string, std::mt19937_64&); + + [[nodiscard]] double getTau() const { return m_tau; } + void setTau(const double tau) { m_tau = tau; } + [[nodiscard]] std::string getName() const { return m_name; }; + + virtual void setParameters(const std::unordered_map&); + + virtual void initialize() = 0; + virtual void draw_tau() = 0; + [[nodiscard]] double tau_wait() const; + + static Distribution* createFromInput(const std::string& input, std::mt19937_64& rng); + +protected: + std::string m_name{"BaseDistribution"}; + double m_tau{1.}; + double m_tau_jump{1.}; + std::mt19937_64& m_rng; +}; + +#endif //RWSIM_TIMESBASE_H diff --git a/times/delta.cpp b/src/times/delta.cpp similarity index 71% rename from times/delta.cpp rename to src/times/delta.cpp index cec764d..9946ff7 100644 --- a/times/delta.cpp +++ b/src/times/delta.cpp @@ -1,10 +1,7 @@ -// -// Created by dominik on 8/12/24. -// #include "delta.h" -DeltaDistribution::DeltaDistribution(const double tau, std::mt19937_64& rng) : Distribution(tau, rng) {} -DeltaDistribution::DeltaDistribution(std::mt19937_64& rng) : Distribution(rng) {} +DeltaDistribution::DeltaDistribution(const double tau, std::mt19937_64& rng) : Distribution(std::string("Delta"), tau, rng) {} +DeltaDistribution::DeltaDistribution(std::mt19937_64& rng) : Distribution(std::string("Delta"), rng) {} void DeltaDistribution::initialize() { m_tau_jump = m_tau; diff --git a/times/delta.h b/src/times/delta.h similarity index 89% rename from times/delta.h rename to src/times/delta.h index 4367db6..463fefc 100644 --- a/times/delta.h +++ b/src/times/delta.h @@ -1,7 +1,3 @@ -// -// Created by dominik on 8/12/24. -// - #ifndef RWSIM_TIMESDELTA_H #define RWSIM_TIMESDELTA_H diff --git a/times/lognormal.cpp b/src/times/lognormal.cpp similarity index 51% rename from times/lognormal.cpp rename to src/times/lognormal.cpp index 8e2ca6f..e6c2b5a 100644 --- a/times/lognormal.cpp +++ b/src/times/lognormal.cpp @@ -1,12 +1,13 @@ -// -// Created by dominik on 8/24/24. -// - #include "lognormal.h" #include -LogNormalDistribution::LogNormalDistribution(const double tau, const double sigma, std::mt19937_64& rng) : Distribution(tau, rng), m_sigma(sigma), m_distribution(std::log(tau), sigma) {} -LogNormalDistribution::LogNormalDistribution(std::mt19937_64& rng) : Distribution(rng) {} +LogNormalDistribution::LogNormalDistribution(const double tau, const double sigma, std::mt19937_64& rng) : Distribution(std::string("LogNormal"), tau, rng), m_sigma(sigma), m_distribution(std::log(tau), sigma) {} +LogNormalDistribution::LogNormalDistribution(std::mt19937_64& rng) : Distribution(std::string("LogNormal"), rng) {} + +void LogNormalDistribution::setParameters(const std::unordered_map ¶meters) { + m_sigma = parameters.at("sigma"); + Distribution::setParameters(parameters); +} void LogNormalDistribution::initialize() { m_distribution = std::lognormal_distribution(std::log(m_tau), m_sigma); diff --git a/times/lognormal.h b/src/times/lognormal.h similarity index 81% rename from times/lognormal.h rename to src/times/lognormal.h index 818c35c..da89503 100644 --- a/times/lognormal.h +++ b/src/times/lognormal.h @@ -1,15 +1,17 @@ - #ifndef LOGNORMAL_H #define LOGNORMAL_H #include "base.h" #include +#include class LogNormalDistribution final : public Distribution { public: LogNormalDistribution(double, double, std::mt19937_64&); explicit LogNormalDistribution(std::mt19937_64 &rng); + void setParameters(const std::unordered_map &) override; + void initialize() override; void draw_tau() override; diff --git a/utils/functions.cpp b/src/utils/functions.cpp similarity index 100% rename from utils/functions.cpp rename to src/utils/functions.cpp diff --git a/utils/functions.h b/src/utils/functions.h similarity index 100% rename from utils/functions.h rename to src/utils/functions.h diff --git a/utils/io.cpp b/src/utils/io.cpp similarity index 61% rename from utils/io.cpp rename to src/utils/io.cpp index a099dfc..b7afbe9 100644 --- a/utils/io.cpp +++ b/src/utils/io.cpp @@ -57,7 +57,7 @@ Arguments parse_args(const int argc, char* argv[]) { continue; } - // all other optional parameter are + // all other arguments are optional parameter auto [option_name, option_value] = get_optional_parameter(it); args.optional[option_name] = option_value; continue; @@ -74,10 +74,15 @@ Arguments parse_args(const int argc, char* argv[]) { continue; } + if (args.distribution_type.empty()) { + args.distribution_type = *it; + continue; + } + throw std::invalid_argument("too many positional arguments"); } - if (args.parameter_file.empty() || args.motion_type.empty()) { + if (args.parameter_file.empty() || args.motion_type.empty() || args.distribution_type.empty()) { throw std::invalid_argument("Missing argument"); } @@ -85,8 +90,6 @@ Arguments parse_args(const int argc, char* argv[]) { } - - std::unordered_map read_parameter(const std::filesystem::path& infile) { if (!std::filesystem::exists(infile)) { std::cerr << "File " << infile << " does not exist" << std::endl; @@ -121,44 +124,72 @@ std::unordered_map read_parameter(const std::filesystem::pa } -void fid_write_out(const std::string& filename, const std::vector& x, const std::vector& y, const double tau, const double t_evo) { - auto size = x.size(); +void save_parameter_to_file( + const std::string& resulttype, + const std::string& motiontype, + const std::string& disttype, + std::unordered_map& parameter, + std::unordered_map& optional + ) { + + std::ostringstream parameter_filename; + parameter_filename << resulttype << "_" << motiontype << "_" << disttype; + + parameter_filename << std::setprecision(6) << std::scientific; + for (const auto& [key, value]: optional) { + parameter_filename << "_" << key << "=" << value; + } + parameter_filename << "_parameter.txt"; - std::ostringstream sfile; - sfile << filename << "_"; - sfile << std::setprecision(6) << std::scientific; - sfile << "tau=" << tau << "_tevo=" << t_evo << ".dat"; { - std::string outfile = sfile.str(); - std::ofstream fid_file(outfile, std::ios::out); - for (unsigned int i = 0; i < size; i++) { - fid_file << x[i] << "\t" << y[i] << "\n"; + // write data to file, columns are secondary axis (echo delay, evolution times) + std::string datafile = parameter_filename.str(); + std::ofstream filestream(datafile, std::ios::out); + + for (const auto& [key, value]: parameter) { + filestream << key << "=" << value << "\n"; } } } -void fid_write_out(const std::string& filename, const std::vector& x, const std::map>& y, const double tau) { - auto size = x.size(); - std::ostringstream sfile; - sfile << filename << "_"; - sfile << std::setprecision(6) << std::scientific; - sfile << "tau=" << tau << ".dat"; +void save_data_to_file( + const std::string& resulttype, + const std::string& motiontype, + const std::string& disttype, + const std::vector& x, + const std::map>& 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"; + { - std::string outfile = sfile.str(); - std::ofstream fid_file(outfile, std::ios::out); - fid_file << "#"; + // 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); + + // first line are values of secondary axis + filestream << "#"; for (const auto& [t_echo_j, _] : y) { - fid_file << t_echo_j << "\t"; + filestream << t_echo_j << "\t"; } - fid_file << std::endl; + filestream << std::endl; + // write values to file + auto size = x.size(); for (unsigned int i = 0; i < size; i++) { - fid_file << x[i]; + filestream << x[i]; for (const auto& [_, fid_j] : y) { - fid_file << "\t" << fid_j[i]; + filestream << "\t" << fid_j[i]; } - fid_file << "\n"; + filestream << "\n"; } } } diff --git a/utils/io.h b/src/utils/io.h similarity index 58% rename from utils/io.h rename to src/utils/io.h index acae839..e26fbf5 100644 --- a/utils/io.h +++ b/src/utils/io.h @@ -1,4 +1,3 @@ - #ifndef RWSIM_IO_H #define RWSIM_IO_H @@ -13,6 +12,7 @@ struct Arguments { bool ste = false; bool spectrum = false; std::string motion_type{}; + std::string distribution_type{}; std::unordered_map optional; }; @@ -21,7 +21,7 @@ std::pair get_optional_parameter(std::vector:: 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); +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&); #endif diff --git a/utils/ranges.cpp b/src/utils/ranges.cpp similarity index 96% rename from utils/ranges.cpp rename to src/utils/ranges.cpp index d10faac..6d16e24 100644 --- a/utils/ranges.cpp +++ b/src/utils/ranges.cpp @@ -1,8 +1,3 @@ -// -// Created by dominik on 8/14/24. -// - - #include #include #include diff --git a/utils/ranges.h b/src/utils/ranges.h similarity index 100% rename from utils/ranges.h rename to src/utils/ranges.h diff --git a/test.py b/test.py deleted file mode 100644 index 4eb567e..0000000 --- a/test.py +++ /dev/null @@ -1,219 +0,0 @@ -import pathlib -import subprocess - -import numpy as np -from scipy.optimize import curve_fit -import matplotlib.pyplot as plt - - - -def run_sims(taus, ste: bool = True, spectrum: bool = False, exec_file: str = './build/rwsim', config_file: str = './config.txt'): - for tau in taus: - arguments = [exec_file, config_file] - if ste: - arguments += ['--ste'] - if spectrum: - arguments += ['--spectrum'] - - arguments += ['BimodalAngle'] - - arguments += ['-TAU', f'{tau}'] - - subprocess.run(arguments) - - - -def dampening(x: np.ndarray, apod: float) -> np.ndarray: - return np.exp(-apod * x) - - -def pulse_attn(freq: np.ndarray, t_pulse: float): - # cf. Schmitt-Rohr/Spieß eq. 2.126; omega_1 * t_p = pi/2 - pi_half_squared = np.pi**2 / 4 - omega = 2 * np.pi * freq - - numerator = np.sin(np.sqrt(pi_half_squared + omega**2 * t_pulse**2 / 2)) - denominator = np.sqrt(pi_half_squared + omega**2 * t_pulse**2 / 4) - - return np.pi * numerator / denominator / 2 - - -def post_process_spectrum(taus, apod, tpulse): - reduction_factor = np.zeros((taus.size, 5)) # hard-coded t_echo :( - - for i, tau in enumerate(taus): - try: - raw_data = np.loadtxt(f'fid_tau={tau:.6e}.dat') - except OSError: - continue - - t = raw_data[:, 0] - timesignal = raw_data[:, 1:] - - timesignal *= dampening(t, apod)[:, None] - timesignal[0, :] /= 2 - - # FT to spectrum - freq = np.fft.fftshift(np.fft.fftfreq(t.size, d=1e-6)) - spec = np.fft.fftshift(np.fft.fft(timesignal, axis=0), axes=0).real - spec *= pulse_attn(freq, t_pulse=tpulse)[:, None] - - reduction_factor[i, :] = 2*timesignal[0, :] - - plt.plot(freq, spec) - plt.show() - - plt.semilogx(taus, reduction_factor, '.') - plt.show() - - -def post_process_ste(taus): - tevo = np.linspace(1e-6, 120e-6, num=121) - - for i, tau in enumerate(taus): - try: - raw_data_cc = np.loadtxt(f'coscos_tau={tau:.6e}.dat') - raw_data_ss = np.loadtxt(f'sinsin_tau={tau:.6e}.dat') - except OSError: - continue - - t_mix = raw_data_cc[:, 0] - - 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:], '.') - - 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:]) - - fig_tau, ax_tau = plt.subplots() - fig_beta, ax_beta= plt.subplots() - fig_finfty, ax_finfty = plt.subplots() - - tau_cc_fit = [] - beta_cc_fit = [] - finfty_cc_fit = [] - - tau_ss_fit = [] - 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)))], - 0.5, - 0.1 - ] - - res = curve_fit(ste, t_mix, 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'Cos-Cos-Fit for {tevo[j]}: tau_c = {tauc:.6e}, beta={beta:.4e}, amplitude = {m0: .4e}, f_infty={finfty:.4e}') - l = ax_cc_raw.semilogx(t_mix, raw_data_cc[:, j+1], 'x', label=f'{tevo[j]}') - ax_cc_raw.semilogx(t_mix, ste(t_mix, *res[0]), linestyle='--', color = l[0].get_color()) - - tau_cc_fit.append(res[0][1]) - beta_cc_fit.append(res[0][2]) - finfty_cc_fit.append(res[0][3]) - - p0 = [ - raw_data_cc[0, 1], - t_mix[np.argmin(np.abs(scaled_ss[:, j]-np.exp(-1)))], - 0.5, - 0.1 - ] - - res = curve_fit(ste, t_mix, raw_data_ss[:, 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_ss_fit.append(res[0][1]) - beta_ss_fit.append(res[0][2]) - finfty_ss_fit.append(res[0][3]) - - p0 = [ - raw_data_cc[0, 1], - t_mix[np.argmin(np.abs(scaled_ss[:, j]-np.exp(-1)))], - 0.5, - 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:], 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-') - - ax_cc_raw.legend() - 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))) - - - - # for i, tau in enumerate(taus): - # - # try: - # raw_data_cc = np.loadtxt(f'coscos_tau={tau:.6e}.dat') - # raw_data_ss = np.loadtxt(f'sinsin_tau={tau:.6e}.dat') - # except OSError: - # continue - # - # t_mix = raw_data_cc[:, 0] - # - # plt.semilogx(t_mix, raw_data_cc[:, 1:]) - # plt.show() - # - # plt.semilogx(t_mix, raw_data_ss[:, 1:]) - # plt.show() - # - # plt.plot(raw_data_cc[0, 1:]) - # plt.plot(raw_data_ss[0, 1:]) - # plt.show() - # - # plt.plot(raw_data_cc[-1, 1:]/raw_data_cc[0, 1:]) - # plt.plot(raw_data_ss[-1, 1:]/raw_data_ss[0, 1:]) - # plt.show() - - -def ste(x, m0, t, beta, finfty): - return m0 * ((1-finfty) * np.exp(-(x/t)**beta) + finfty) - - -if __name__ == '__main__': - tau_values = np.geomspace(1e-2, 1e2, num=2) # if num=1, tau is first value - - # parameter for spectrum simulations - lb = 2e3 - pulse_length = 2e-6 - - run_sims(tau_values) - post_process_ste(tau_values) - # post_process_spectrum(tau_values, lb, pulse_length) \ No newline at end of file diff --git a/times/base.cpp b/times/base.cpp deleted file mode 100644 index b14e967..0000000 --- a/times/base.cpp +++ /dev/null @@ -1,13 +0,0 @@ -// -// Created by dominik on 8/12/24. -// -#include "base.h" - -Distribution::Distribution(const double tau, std::mt19937_64 &rng) : m_tau(tau), m_tau_jump(tau), m_rng(rng) {} - -Distribution::Distribution(std::mt19937_64 &rng) : m_rng(rng) {} - -double Distribution::tau_wait() const { - return std::exponential_distribution(1./m_tau_jump)(m_rng); -} - diff --git a/times/base.h b/times/base.h deleted file mode 100644 index bdd1e0b..0000000 --- a/times/base.h +++ /dev/null @@ -1,30 +0,0 @@ -// -// Created by dominik on 8/12/24. -// - -#ifndef RWSIM_TIMESBASE_H -#define RWSIM_TIMESBASE_H - -#include - -class Distribution { -public: - virtual ~Distribution() = default; - - Distribution(double, std::mt19937_64&); - explicit Distribution(std::mt19937_64&); - - [[nodiscard]] double getTau() const { return m_tau; } - void setTau(const double tau) { m_tau = tau;} - - virtual void initialize() = 0; - virtual void draw_tau() = 0; - [[nodiscard]] double tau_wait() const; - -protected: - double m_tau{1.}; - double m_tau_jump{1.}; - std::mt19937_64& m_rng; -}; - -#endif //RWSIM_TIMESBASE_H