diff --git a/README.md b/README.md index cfd88b2..91bf211 100644 --- a/README.md +++ b/README.md @@ -1,32 +1,37 @@ -# Download +# NMR Random Walk Simulation -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` +Simulate solid-state NMR spectra and stimulated echo (STE) decays using random walk models. # Build -# Running - -## Command line - -To run a random walk simulation Execute `rwsim` in the command line with +Requires CMake 3.18+, a C++17 compiler, and OpenMP. ```bash - ./rwsim /path/to/config.txt MotionModel DistributionModel --ste --spectrum -ARG1 1 -ARG2 2 +cmake -B build +cmake --build build ``` -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. +The executable is `build/src/rwsim`. +# Running -## Configuration file `config.txt` +```bash +./build/src/rwsim /path/to/config.txt MotionModel DistributionModel --ste --spectrum -ARG1 1 -ARG2 2 +``` -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. +Three positional arguments are required: the path to a configuration file, the motion model name, and the distribution model name. +Set `--ste` and/or `--spectrum` to run stimulated echo or spectrum simulations. +Any parameter from the configuration file can be overridden on the command line, e.g. `-tau 1e-3`. -### General parameter +Use `-seed 42` for reproducible results. Without a seed, a random one is used. -This list of general parameter are necessary for all simulations: +**Output files are overwritten if a simulation with the same parameters is run again.** + +## Configuration file + +Change delta, eta, mixing times, echo times, etc. in this file. If a parameter appears multiple times, only the last value is used. + +### General parameters | Parameter | Description | |------------|----------------------------| @@ -36,48 +41,88 @@ This list of general parameter are necessary for all simulations: ### 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 | +| Model | CLI name | Parameters | +|-------|----------|------------| +| Delta distribution (same tau for all walkers) | `Delta` | `tau` (jump time in s) | +| Log-normal distribution | `LogNormal` | `tau` (peak time in s), `sigma` (width) | ### 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`. +| Model | CLI name | Parameters | +|-------|----------|------------| +| Isotropic random jump | `RandomJump` | — | +| Isotropic jump by fixed angle | `IsotropicAngle` | `angle` (degrees) | +| Bimodal angle distribution | `BimodalAngle` | `angle1`, `angle2` (degrees), `probability1` (0–1) | +| Tetrahedral 4-site jump | `FourSiteTetrahedral` | — | +| Octahedral 6-site jump (C3) | `SixSiteOctahedralC3` | — | +| N-site jump on cone | `NSiteConeJump` | `num_sites`, `chi` (cone angle, degrees) | +| Random jump on cone | `RandomJumpOnCone` | `angle` (cone angle) | +| Wobbling in cone | `ConeWobble` | `angle` (cone angle) | -+ 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` +### Spectrum parameters - | Parameter | Description | - |-----------|----------------------| - | angle | Jump angle in degree | +| Parameter | Description | +|--------------|---------------------------------| +| dwell_time | Acquisition dwell time in s | +| num_acq | Number of acquisition points | +| techo_start | First echo time in s | +| techo_stop | Last echo time in s | +| techo_steps | Number of echo times | -+ Parameters for isotropic jump with bimodal angle distribution `BimodalAngle` +### STE parameters - | Parameter | Description | - |--------------|--------------------------------------| - | angle1 | First jump angle in degree | - | angle2 | Second jump angle in degree | - | probability1 | Probality that jump has angle1 (0-1) | +| Parameter | Description | +|-------------|---------------------------------| +| tevo_start | First evolution time in s | +| tevo_stop | Last evolution time in s | +| tevo_steps | Number of evolution times | +| tmix_start | First mixing time in s | +| tmix_stop | Last mixing time in s | +| tmix_steps | Number of mixing times (log-spaced) | +| tpulse4 | Delay of 4th pulse in s | +# Architecture +The simulation is built around three abstractions: +## Dynamics +`Dynamics` (`src/dynamics/base.h`) produces trajectory steps — pairs of `{dt, omega}` representing a waiting time and the NMR frequency after a jump. This is the core interface for defining physical models. +**`DecoupledDynamics`** wraps a separate motion model and time distribution, calling them independently. This covers all existing models. -# Running +To implement a **coupled model** where motion and time are interdependent (e.g. jump angle affects waiting time, or position-dependent dynamics), implement the `Dynamics` interface directly: -**Because filenames are always the same, previous simulations and results are overwritten!** +```cpp +class MyModel : public Dynamics { + void initialize(std::mt19937_64& rng) override { /* set initial state */ } + Step next(std::mt19937_64& rng) override { + // draw angle, compute rate from angle, draw wait time from rate + return {dt, omega}; + } + // ... clone(), setParameters(), etc. +}; +``` -1. Execute `build.sh` (in terminal with `./build.sh`). It compiles the source code and executes test.py +## Experiment + +`Experiment` (`src/experiments/base.h`) defines how trajectory data is turned into observables: + +- `setup(params)` — configure time axes, allocate accumulators +- `accumulate(trajectory)` — process one walker's trajectory +- `save(directory)` — write results to files +- `clone()` / `merge()` — for thread-parallel accumulation + +Existing experiments: `SpectrumExperiment`, `STEExperiment`. + +## Simulation runner + +`run_simulation()` (`src/sims.h`) connects a `Dynamics` and an `Experiment`: + +1. Sets parameters on both +2. Generates trajectories in parallel (OpenMP, one clone per thread) +3. Each thread accumulates into its own experiment clone +4. Merges thread-local results +5. Saves output + +The walker loop is parallelized with OpenMP using static scheduling for deterministic reproducibility with a given seed and thread count. diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index bbfa3b5..5f861d5 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -3,15 +3,16 @@ add_subdirectory(times) add_subdirectory(motions) add_subdirectory(utils) add_subdirectory(experiments) +add_subdirectory(dynamics) add_library(simulation STATIC sims.cpp sims.h) -target_link_libraries(simulation PRIVATE utils experiments OpenMP::OpenMP_CXX) +target_link_libraries(simulation PRIVATE utils experiments dynamics OpenMP::OpenMP_CXX) add_executable( rwsim main.cpp ) -target_link_libraries(rwsim PUBLIC RWMotion RWTime utils experiments simulation) +target_link_libraries(rwsim PUBLIC RWMotion RWTime utils experiments dynamics simulation) target_compile_options(rwsim PUBLIC -Werror -Wall -Wextra -Wconversion -O2) diff --git a/src/dynamics/CMakeLists.txt b/src/dynamics/CMakeLists.txt new file mode 100644 index 0000000..ec675d2 --- /dev/null +++ b/src/dynamics/CMakeLists.txt @@ -0,0 +1,6 @@ +add_library( + dynamics STATIC + base.h + decoupled.cpp + decoupled.h +) diff --git a/src/dynamics/base.h b/src/dynamics/base.h new file mode 100644 index 0000000..a71e62c --- /dev/null +++ b/src/dynamics/base.h @@ -0,0 +1,26 @@ +#ifndef RWSIM_DYNAMICSBASE_H +#define RWSIM_DYNAMICSBASE_H + +#include +#include +#include +#include + +struct Step { + double dt; + double omega; +}; + +class Dynamics { +public: + virtual ~Dynamics() = default; + + virtual void setParameters(const std::unordered_map& parameters) = 0; + virtual void initialize(std::mt19937_64& rng) = 0; + virtual Step next(std::mt19937_64& rng) = 0; + [[nodiscard]] virtual double getInitOmega() const = 0; + [[nodiscard]] virtual std::unique_ptr clone() const = 0; + [[nodiscard]] virtual std::string toString() const = 0; +}; + +#endif //RWSIM_DYNAMICSBASE_H diff --git a/src/dynamics/decoupled.cpp b/src/dynamics/decoupled.cpp new file mode 100644 index 0000000..681c4cb --- /dev/null +++ b/src/dynamics/decoupled.cpp @@ -0,0 +1,35 @@ +#include "decoupled.h" + +DecoupledDynamics::DecoupledDynamics( + std::unique_ptr motion, + std::unique_ptr dist) + : m_motion(std::move(motion)), m_dist(std::move(dist)) {} + +void DecoupledDynamics::setParameters(const std::unordered_map& parameters) { + m_motion->setParameters(parameters); + m_dist->setParameters(parameters); +} + +void DecoupledDynamics::initialize(std::mt19937_64& rng) { + m_motion->initialize(rng); + m_dist->initialize(rng); +} + +Step DecoupledDynamics::next(std::mt19937_64& rng) { + double dt = m_dist->tau_wait(rng); + double omega = m_motion->jump(rng); + return {dt, omega}; +} + +double DecoupledDynamics::getInitOmega() const { + return m_motion->getInitOmega(); +} + +std::unique_ptr DecoupledDynamics::clone() const { + return std::make_unique( + m_motion->clone(), m_dist->clone()); +} + +std::string DecoupledDynamics::toString() const { + return m_motion->toString() + "/" + m_dist->toString(); +} diff --git a/src/dynamics/decoupled.h b/src/dynamics/decoupled.h new file mode 100644 index 0000000..6be181d --- /dev/null +++ b/src/dynamics/decoupled.h @@ -0,0 +1,27 @@ +#ifndef RWSIM_DECOUPLED_H +#define RWSIM_DECOUPLED_H + +#include "base.h" +#include "../motions/base.h" +#include "../times/base.h" + +#include + +class DecoupledDynamics final : public Dynamics { +public: + DecoupledDynamics(std::unique_ptr motion, + std::unique_ptr dist); + + void setParameters(const std::unordered_map& parameters) override; + void initialize(std::mt19937_64& rng) override; + Step next(std::mt19937_64& rng) override; + [[nodiscard]] double getInitOmega() const override; + [[nodiscard]] std::unique_ptr clone() const override; + [[nodiscard]] std::string toString() const override; + +private: + std::unique_ptr m_motion; + std::unique_ptr m_dist; +}; + +#endif //RWSIM_DECOUPLED_H diff --git a/src/experiments/base.h b/src/experiments/base.h index 2688545..d13c34e 100644 --- a/src/experiments/base.h +++ b/src/experiments/base.h @@ -1,9 +1,6 @@ #ifndef RWSIM_EXPERIMENTBASE_H #define RWSIM_EXPERIMENTBASE_H -#include "../motions/base.h" -#include "../times/base.h" - #include #include #include @@ -23,7 +20,7 @@ public: const std::unordered_map& optional) = 0; [[nodiscard]] virtual double tmax() const = 0; virtual void accumulate(const Trajectory& traj, double init_omega, int num_walkers) = 0; - virtual void save(const motions::BaseMotion& motion, const times::BaseDistribution& dist) = 0; + virtual void save(const std::string& directory) = 0; [[nodiscard]] virtual std::unique_ptr clone() const = 0; virtual void merge(const Experiment& other) = 0; }; diff --git a/src/experiments/spectrum.cpp b/src/experiments/spectrum.cpp index 8477ce6..6439db4 100644 --- a/src/experiments/spectrum.cpp +++ b/src/experiments/spectrum.cpp @@ -3,6 +3,7 @@ #include "../utils/ranges.h" #include "../utils/io.h" + #include #include @@ -44,10 +45,9 @@ void SpectrumExperiment::accumulate(const Trajectory& traj, double, int num_walk } } -void SpectrumExperiment::save(const motions::BaseMotion& motion, const times::BaseDistribution& dist) { - const auto path = make_directory(motion, dist); - save_parameter_to_file(std::string("timesignal"), path, m_parameter, m_optional); - save_data_to_file(std::string("timesignal"), path, m_t_fid, m_fid_dict, m_optional); +void SpectrumExperiment::save(const std::string& directory) { + save_parameter_to_file(std::string("timesignal"), directory, m_parameter, m_optional); + save_data_to_file(std::string("timesignal"), directory, m_t_fid, m_fid_dict, m_optional); } std::unique_ptr SpectrumExperiment::clone() const { diff --git a/src/experiments/spectrum.h b/src/experiments/spectrum.h index 2ac40f2..86155c9 100644 --- a/src/experiments/spectrum.h +++ b/src/experiments/spectrum.h @@ -13,7 +13,7 @@ public: const std::unordered_map& optional) override; [[nodiscard]] double tmax() const override; void accumulate(const Trajectory& traj, double init_omega, int num_walkers) override; - void save(const motions::BaseMotion& motion, const times::BaseDistribution& dist) override; + void save(const std::string& directory) override; [[nodiscard]] std::unique_ptr clone() const override; void merge(const Experiment& other) override; diff --git a/src/experiments/ste.cpp b/src/experiments/ste.cpp index d1224df..7808ac5 100644 --- a/src/experiments/ste.cpp +++ b/src/experiments/ste.cpp @@ -70,12 +70,11 @@ void STEExperiment::accumulate(const Trajectory& traj, double init_omega, int nu } } -void STEExperiment::save(const motions::BaseMotion& motion, const times::BaseDistribution& dist) { - const auto folders = make_directory(motion, dist); - save_parameter_to_file(std::string("ste"), folders, m_parameter, m_optional); - save_data_to_file(std::string("coscos"), folders, m_mixing_times, m_cc_dict, m_optional); - save_data_to_file(std::string("sinsin"), folders, m_mixing_times, m_ss_dict, m_optional); - save_data_to_file(std::string("f2"), folders, m_mixing_times, m_f2, m_optional); +void STEExperiment::save(const std::string& directory) { + save_parameter_to_file(std::string("ste"), directory, m_parameter, m_optional); + save_data_to_file(std::string("coscos"), directory, m_mixing_times, m_cc_dict, m_optional); + save_data_to_file(std::string("sinsin"), directory, m_mixing_times, m_ss_dict, m_optional); + save_data_to_file(std::string("f2"), directory, m_mixing_times, m_f2, m_optional); } std::unique_ptr STEExperiment::clone() const { diff --git a/src/experiments/ste.h b/src/experiments/ste.h index a276d77..09766cc 100644 --- a/src/experiments/ste.h +++ b/src/experiments/ste.h @@ -13,7 +13,7 @@ public: const std::unordered_map& optional) override; [[nodiscard]] double tmax() const override; void accumulate(const Trajectory& traj, double init_omega, int num_walkers) override; - void save(const motions::BaseMotion& motion, const times::BaseDistribution& dist) override; + void save(const std::string& directory) override; [[nodiscard]] std::unique_ptr clone() const override; void merge(const Experiment& other) override; diff --git a/src/main.cpp b/src/main.cpp index 02861a7..ad53db3 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,5 +1,6 @@ #include "sims.h" +#include "dynamics/decoupled.h" #include "experiments/spectrum.h" #include "experiments/ste.h" #include "utils/io.h" @@ -43,16 +44,17 @@ int main (const int argc, char *argv[]) { rng.seed(rd()); } - auto motion = motions::BaseMotion::createFromInput(args.motion_type); - auto dist = times::BaseDistribution::createFromInput(args.distribution_type); + DecoupledDynamics dynamics( + motions::BaseMotion::createFromInput(args.motion_type), + times::BaseDistribution::createFromInput(args.distribution_type)); if (args.spectrum) { SpectrumExperiment experiment; - run_simulation(experiment, parameter, args.optional, *motion, *dist, rng); + run_simulation(experiment, parameter, args.optional, dynamics, rng); } if (args.ste) { STEExperiment experiment; - run_simulation(experiment, parameter, args.optional, *motion, *dist, rng); + run_simulation(experiment, parameter, args.optional, dynamics, rng); } return 0; diff --git a/src/sims.cpp b/src/sims.cpp index f85773b..1fb3b86 100644 --- a/src/sims.cpp +++ b/src/sims.cpp @@ -1,5 +1,6 @@ #include "sims.h" #include "utils/functions.h" +#include "utils/io.h" #include #include @@ -10,14 +11,12 @@ void run_simulation( Experiment& experiment, std::unordered_map& parameter, std::unordered_map& optional, - motions::BaseMotion& motion, - times::BaseDistribution& dist, + Dynamics& dynamics, std::mt19937_64& rng ) { const int num_walker = static_cast(parameter["num_walker"]); - dist.setParameters(parameter); - motion.setParameters(parameter); + dynamics.setParameters(parameter); experiment.setup(parameter, optional); const auto start = printStart(optional); @@ -31,13 +30,11 @@ void run_simulation( thread_rngs.emplace_back(rng()); } - // Create per-thread clones of motion, distribution, and experiment - std::vector> thread_motions; - std::vector> thread_dists; + // Create per-thread clones of dynamics and experiment + std::vector> thread_dynamics; std::vector> thread_experiments; for (int i = 0; i < num_threads; i++) { - thread_motions.push_back(motion.clone()); - thread_dists.push_back(dist.clone()); + thread_dynamics.push_back(dynamics.clone()); thread_experiments.push_back(experiment.clone()); } @@ -48,14 +45,13 @@ void run_simulation( { const int tid = omp_get_thread_num(); auto& local_rng = thread_rngs[tid]; - auto& local_motion = *thread_motions[tid]; - auto& local_dist = *thread_dists[tid]; + auto& local_dynamics = *thread_dynamics[tid]; auto& local_experiment = *thread_experiments[tid]; #pragma omp for schedule(static) for (int mol_i = 0; mol_i < num_walker; mol_i++) { - auto traj = make_trajectory(local_motion, local_dist, experiment.tmax(), local_rng); - local_experiment.accumulate(traj, local_motion.getInitOmega(), num_walker); + auto traj = make_trajectory(local_dynamics, experiment.tmax(), local_rng); + local_experiment.accumulate(traj, local_dynamics.getInitOmega(), num_walker); if (tid == 0) { #pragma omp atomic @@ -73,24 +69,23 @@ void run_simulation( experiment.merge(*thread_experiments[i]); } - experiment.save(motion, dist); + const auto directory = make_directory(dynamics.toString()); + experiment.save(directory); printEnd(start); } Trajectory make_trajectory( - motions::BaseMotion& motion, - times::BaseDistribution& dist, + Dynamics& dynamics, const double t_max, std::mt19937_64& rng ) { double t_passed = 0; double phase = 0; - motion.initialize(rng); - dist.initialize(rng); + dynamics.initialize(rng); - double omega = motion.getInitOmega(); + double omega = dynamics.getInitOmega(); Trajectory traj; traj.time.emplace_back(t_passed); @@ -98,10 +93,10 @@ Trajectory make_trajectory( traj.omega.emplace_back(omega); while (t_passed < t_max) { - const double t = dist.tau_wait(rng); - t_passed += t; - omega = motion.jump(rng); - phase += omega * t; + auto [dt, new_omega] = dynamics.next(rng); + phase += omega * dt; + t_passed += dt; + omega = new_omega; traj.time.emplace_back(t_passed); traj.phase.emplace_back(phase); diff --git a/src/sims.h b/src/sims.h index e64d319..eb96e7e 100644 --- a/src/sims.h +++ b/src/sims.h @@ -1,9 +1,8 @@ #ifndef RWSIM_SIMS_H #define RWSIM_SIMS_H +#include "dynamics/base.h" #include "experiments/base.h" -#include "motions/base.h" -#include "times/base.h" #include #include @@ -14,11 +13,10 @@ void run_simulation( Experiment& experiment, std::unordered_map& parameter, std::unordered_map& optional, - motions::BaseMotion& motion, - times::BaseDistribution& dist, + Dynamics& dynamics, std::mt19937_64& rng); -Trajectory make_trajectory(motions::BaseMotion& motion, times::BaseDistribution& dist, double t_max, std::mt19937_64& rng); +Trajectory make_trajectory(Dynamics& dynamics, double t_max, std::mt19937_64& rng); 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 d0fa0fa..3a67a1f 100644 --- a/src/utils/io.cpp +++ b/src/utils/io.cpp @@ -11,11 +11,6 @@ #include #include -#include "../motions/base.h" -#include "../times/base.h" - - -class Motion; std::pair get_optional_parameter(std::vector::const_iterator &it) { std::string stripped_arg; @@ -128,19 +123,12 @@ std::unordered_map read_parameter(const std::filesystem::pa } -std::string make_directory( - const motions::BaseMotion& motion, - const times::BaseDistribution& distribution - ) { - std::ostringstream path_name; - path_name << motion.toString() << "/" << distribution.toString(); - - if (!std::filesystem::create_directories(path_name.str())) { - std::cout << "Created directory " << path_name.str() << std::endl; +std::string make_directory(const std::string& path_name) { + if (!std::filesystem::create_directories(path_name)) { + std::cout << "Created directory " << path_name << std::endl; } - - return path_name.str(); + return path_name; } diff --git a/src/utils/io.h b/src/utils/io.h index e2493b1..50ab2e4 100644 --- a/src/utils/io.h +++ b/src/utils/io.h @@ -7,9 +7,6 @@ #include #include -#include "../motions/base.h" -#include "../times/base.h" - struct Arguments { std::string parameter_file{}; bool ste = false; @@ -23,7 +20,7 @@ Arguments parse_args(int argc, char* argv[]); std::pair get_optional_parameter(std::vector::const_iterator &it); std::unordered_map read_parameter(const std::filesystem::path&); -std::string make_directory(const motions::BaseMotion&, const times::BaseDistribution&); +std::string make_directory(const std::string& path_name); void save_parameter_to_file(const std::string&, const std::string&, const std::unordered_map&, const std::unordered_map&); void save_data_to_file(const std::string&, const std::string&, const std::vector&, const std::map>&, const std::unordered_map&);