From 5acbaaa5f8f026930641eb36d3925db39da1ae78 Mon Sep 17 00:00:00 2001 From: Dominik Demuth Date: Sun, 8 Mar 2026 15:33:56 +0100 Subject: [PATCH] Comments --- src/dynamics/base.h | 20 ++++++++++++++++++++ src/dynamics/decoupled.h | 4 ++++ src/experiments/base.h | 23 ++++++++++++++++++++--- src/experiments/spectrum.h | 16 ++++++++++++---- src/experiments/ste.h | 22 ++++++++++++++++------ src/motions/base.cpp | 2 ++ src/motions/base.h | 26 +++++++++++++++++++++++--- src/motions/conemotion.h | 9 +++++++-- src/motions/coordinates.h | 8 +++++++- src/motions/diffusivemotion.h | 7 ++++++- src/sims.h | 7 +++++++ src/times/base.h | 22 ++++++++++++++++++---- src/utils/functions.cpp | 7 ++++--- src/utils/registry.h | 11 +++++++++++ 14 files changed, 157 insertions(+), 27 deletions(-) diff --git a/src/dynamics/base.h b/src/dynamics/base.h index 4b38e11..e0ec08a 100644 --- a/src/dynamics/base.h +++ b/src/dynamics/base.h @@ -6,23 +6,43 @@ #include #include +/// A single step in the dynamics: waiting time dt and new frequency omega. struct Step { double dt; double omega; }; +/// Interface for trajectory dynamics (motion + timing). +/// +/// Produces a sequence of (dt, omega) steps that define a walker's trajectory. +/// DecoupledDynamics wraps independent motion and time models. +/// Coupled models (e.g., position-dependent jump rates) implement this directly. +/// +/// Use 2 positional CLI args for coupled: ./rwsim config.txt MyDynamics +/// Use 3 positional CLI args for decoupled: ./rwsim config.txt Motion Dist class Dynamics { public: virtual ~Dynamics() = default; virtual void setParameters(const std::unordered_map ¶meters) = 0; + + /// Set up initial state for a new walker (e.g., draw random orientation). virtual void initialize(std::mt19937_64 &rng) = 0; + + /// Generate the next step: waiting time and new NMR frequency. virtual Step next(std::mt19937_64 &rng) = 0; + + /// Return the NMR frequency at the walker's initial position. [[nodiscard]] virtual double getInitOmega() const = 0; + + /// Deep copy for per-thread cloning in parallel simulation. [[nodiscard]] virtual std::unique_ptr clone() const = 0; + + /// String identifier used for output directory naming. [[nodiscard]] virtual std::string toString() const = 0; + /// Factory: create a Dynamics instance by registered name. static std::unique_ptr createFromInput(const std::string &input); }; diff --git a/src/dynamics/decoupled.h b/src/dynamics/decoupled.h index 81d4d7a..44e145a 100644 --- a/src/dynamics/decoupled.h +++ b/src/dynamics/decoupled.h @@ -7,6 +7,10 @@ #include +/// Dynamics where motion and waiting times are independent. +/// Wraps a BaseMotion (geometry) and BaseDistribution (timing) and delegates +/// to each independently. This is the default mode when two positional CLI +/// args are given (motion + distribution). class DecoupledDynamics final : public Dynamics { public: DecoupledDynamics(std::unique_ptr motion, diff --git a/src/experiments/base.h b/src/experiments/base.h index 118a734..98cfc01 100644 --- a/src/experiments/base.h +++ b/src/experiments/base.h @@ -6,24 +6,41 @@ #include #include +/// Time series of a single walker: NMR frequency, accumulated phase, and time. struct Trajectory { - std::vector time; - std::vector phase; - std::vector omega; + std::vector time; ///< Cumulative time at each step + std::vector phase; ///< Accumulated NMR phase (integral of omega*dt) + std::vector omega; ///< Instantaneous NMR frequency at each step }; +/// Interface for NMR experiments that process walker trajectories. +/// +/// Lifecycle: setup() -> [accumulate() per walker] -> merge() -> save() +/// Each thread gets a clone(); results are merged after the parallel loop. class Experiment { public: virtual ~Experiment() = default; + /// Configure time axes and allocate output buffers from simulation parameters. virtual void setup(const std::unordered_map ¶meter, const std::unordered_map &optional) = 0; + + /// Maximum trajectory time needed by this experiment. [[nodiscard]] virtual double tmax() const = 0; + + /// Process one walker's trajectory, adding its contribution to the result. + /// Results are normalized by num_walkers during accumulation. virtual void accumulate(const Trajectory &traj, double init_omega, int num_walkers) = 0; + + /// Write results to files in the given directory. virtual void save(const std::string &directory) = 0; + + /// Deep copy for per-thread cloning in parallel simulation. [[nodiscard]] virtual std::unique_ptr clone() const = 0; + + /// Combine results from another (per-thread) experiment instance. virtual void merge(const Experiment &other) = 0; }; diff --git a/src/experiments/spectrum.h b/src/experiments/spectrum.h index 22dfd75..d5c0e89 100644 --- a/src/experiments/spectrum.h +++ b/src/experiments/spectrum.h @@ -7,6 +7,14 @@ #include #include +/// Solid-echo NMR spectrum experiment. +/// +/// Computes the free induction decay (FID) after a solid echo at various echo +/// times. For each echo time t_echo, the FID is sampled at t_fid points after +/// refocusing at 2*t_echo. The signal is cos(phase_acq - 2*phase_echo), +/// averaged over all walkers. +/// +/// Output: timesignal files with FID data for each echo time. class SpectrumExperiment final : public Experiment { public: void setup(const std::unordered_map ¶meter, @@ -19,10 +27,10 @@ public: void merge(const Experiment &other) override; private: - int m_num_acq{0}; - std::vector m_t_fid; - std::vector m_echo_times; - std::map> m_fid_dict; + int m_num_acq{0}; ///< Number of acquisition points per FID + std::vector m_t_fid; ///< FID time axis (dwell_time spacing) + std::vector m_echo_times; ///< Echo delay times + std::map> m_fid_dict; ///< FID data per echo time std::unordered_map m_parameter; std::unordered_map m_optional; double m_tmax{0}; diff --git a/src/experiments/ste.h b/src/experiments/ste.h index a4f831a..92147bd 100644 --- a/src/experiments/ste.h +++ b/src/experiments/ste.h @@ -7,6 +7,16 @@ #include #include +/// Stimulated echo (STE) experiment. +/// +/// Computes the stimulated echo decay as a function of mixing time for various +/// evolution times. The STE probes slow molecular reorientation by comparing +/// frequencies before and after a mixing period. +/// +/// Output files: +/// coscos — cos-cos correlation function F_cc(t_evo, t_mix) +/// sinsin — sin-sin correlation function F_ss(t_evo, t_mix) +/// f2 — two-time correlation function class STEExperiment final : public Experiment { public: void setup(const std::unordered_map ¶meter, @@ -20,12 +30,12 @@ public: private: int m_num_mix_times{0}; - double m_tpulse4{0}; - std::vector m_evolution_times; - std::vector m_mixing_times; - std::map> m_cc_dict; - std::map> m_ss_dict; - std::vector m_f2; + double m_tpulse4{0}; ///< Duration of the 4th pulse + std::vector m_evolution_times; ///< t_evo values (linear) + std::vector m_mixing_times; ///< t_mix values (logarithmic) + std::map> m_cc_dict; ///< cos-cos per t_evo + std::map> m_ss_dict; ///< sin-sin per t_evo + std::vector m_f2; ///< Two-time correlator std::unordered_map m_parameter; std::unordered_map m_optional; double m_tmax{0}; diff --git a/src/motions/base.cpp b/src/motions/base.cpp index ff81e8a..c5d9750 100644 --- a/src/motions/base.cpp +++ b/src/motions/base.cpp @@ -10,6 +10,8 @@ BaseMotion::BaseMotion(std::string name, const double delta, const double eta) BaseMotion::BaseMotion(std::string name) : m_name(std::move(name)) {} +// Quadrupolar NMR frequency: omega_q(theta, phi) = +// pi * delta * (3*cos^2(theta) - 1 - eta*sin^2(theta)*cos(2*phi)) double BaseMotion::omega_q(const double cos_theta, const double phi) const { const double cos_theta_square = cos_theta * cos_theta; const double sin_theta_square = 1. - cos_theta_square; diff --git a/src/motions/base.h b/src/motions/base.h index 91500e1..ce5212a 100644 --- a/src/motions/base.h +++ b/src/motions/base.h @@ -8,6 +8,17 @@ #include namespace motions { + +/// Base class for NMR motion models. +/// +/// Each model defines how a molecular orientation evolves over discrete jumps. +/// The quadrupolar frequency omega_q depends on orientation relative to the +/// static magnetic field via the second-order Legendre polynomial: +/// omega_q = pi * delta * (3*cos^2(theta) - 1 - eta*sin^2(theta)*cos(2*phi)) +/// +/// Intermediate base classes: +/// ConeMotion — models with a fixed symmetry axis (wobble, random jump on cone) +/// DiffusiveMotion — models with continuous rotational diffusion (small angle, bimodal) class BaseMotion { public: virtual ~BaseMotion() = default; @@ -15,12 +26,20 @@ public: BaseMotion(std::string, double, double); explicit BaseMotion(std::string); + /// Draw a uniform random orientation on the unit sphere. coordinates::SphericalPos draw_position(std::mt19937_64 &rng); + + /// Compute quadrupolar frequency from orientation angles. [[nodiscard]] double omega_q(double, double) const; [[nodiscard]] double omega_q(const coordinates::SphericalPos &) const; + /// Set up initial state for a new walker (called once per walker). virtual void initialize(std::mt19937_64 &rng) = 0; + + /// Perform one jump, returning the new omega_q. virtual double jump(std::mt19937_64 &rng) = 0; + + /// Deep copy for per-thread cloning. [[nodiscard]] virtual std::unique_ptr clone() const = 0; virtual void setParameters(const std::unordered_map &); @@ -35,14 +54,15 @@ public: [[nodiscard]] virtual std::string toString() const = 0; + /// Factory: create a motion model by registered name. static std::unique_ptr createFromInput(const std::string &input); protected: std::string m_name{"BaseMotion"}; - double m_delta{1.}; - double m_eta{0.}; + double m_delta{1.}; ///< Quadrupolar coupling constant (Hz) + double m_eta{0.}; ///< Asymmetry parameter [0, 1] std::uniform_real_distribution<> m_uni_dist{0., 1.}; - double m_initial_omega{0.}; + double m_initial_omega{0.}; ///< Frequency at initial position }; std::ostream &operator<<(std::ostream &os, const BaseMotion &m); diff --git a/src/motions/conemotion.h b/src/motions/conemotion.h index 76a8e1b..890adfd 100644 --- a/src/motions/conemotion.h +++ b/src/motions/conemotion.h @@ -5,10 +5,15 @@ #include "coordinates.h" namespace motions { + +/// Intermediate base for motions with a fixed symmetry axis and cone angle. +/// Provides shared axis initialization and angle parameter handling. +/// Used by WobbleCone and RandomJumpOnCone. class ConeMotion : public BaseMotion { public: using BaseMotion::BaseMotion; + /// Draws a random symmetry axis on the unit sphere. void initialize(std::mt19937_64 &rng) override; void setParameters(const std::unordered_map &) override; @@ -16,8 +21,8 @@ public: getParameters() const override; protected: - double m_angle{0}; - coordinates::SphericalPos m_axis{1, 0}; + double m_angle{0}; ///< Half-opening angle of the cone (radians) + coordinates::SphericalPos m_axis{1, 0}; ///< Symmetry axis orientation }; } // namespace motions diff --git a/src/motions/coordinates.h b/src/motions/coordinates.h index 1b1dc9c..6dbfc32 100644 --- a/src/motions/coordinates.h +++ b/src/motions/coordinates.h @@ -2,7 +2,10 @@ #ifndef COORDINATES_H #define COORDINATES_H +/// Coordinate types and rotations for molecular orientations on the unit sphere. namespace coordinates { + +/// Spherical position stored as (cos_theta, phi) to avoid repeated cos/acos. struct SphericalPos { double cos_theta; double phi; @@ -14,7 +17,10 @@ struct CartesianPos { double z; }; -SphericalPos rotate(const SphericalPos &, double, double); +/// Rotate a spherical position by polar angle alpha around an axis defined +/// by azimuthal angle beta (rotation about the original position's z-axis). +SphericalPos rotate(const SphericalPos &, double alpha, double beta); + CartesianPos spherical_to_xyz(const SphericalPos &); SphericalPos xyz_to_spherical(const CartesianPos &); } // namespace coordinates diff --git a/src/motions/diffusivemotion.h b/src/motions/diffusivemotion.h index 0ec37c2..b6244f8 100644 --- a/src/motions/diffusivemotion.h +++ b/src/motions/diffusivemotion.h @@ -5,14 +5,19 @@ #include "coordinates.h" namespace motions { + +/// Intermediate base for motions with continuous rotational diffusion. +/// Tracks the previous position and provides shared initialization. +/// Used by SmallAngle and BimodalAngle. class DiffusiveMotion : public BaseMotion { public: using BaseMotion::BaseMotion; + /// Draws a random initial position and sets the initial frequency. void initialize(std::mt19937_64 &rng) override; protected: - coordinates::SphericalPos m_prev_pos{0., 0.}; + coordinates::SphericalPos m_prev_pos{0., 0.}; ///< Current orientation }; } // namespace motions diff --git a/src/sims.h b/src/sims.h index 41e54c9..14e8435 100644 --- a/src/sims.h +++ b/src/sims.h @@ -10,11 +10,18 @@ #include #include +/// Run a parallel random walk simulation over multiple experiments. +/// +/// Generates one trajectory per walker (using the longest tmax across all +/// experiments) and accumulates results into each experiment. Parallelized +/// with OpenMP using per-thread clones of dynamics and experiments. void run_simulation(std::vector experiments, const std::unordered_map ¶meter, const std::unordered_map &optional, Dynamics &dynamics, std::mt19937_64 &rng); +/// Generate a trajectory by repeatedly calling dynamics.next() until t_max. +/// Reuses the trajectory's allocated memory across calls (clear preserves capacity). void make_trajectory(Trajectory &traj, Dynamics &dynamics, double t_max, std::mt19937_64 &rng); diff --git a/src/times/base.h b/src/times/base.h index bbe7aa0..531d6ce 100644 --- a/src/times/base.h +++ b/src/times/base.h @@ -6,6 +6,16 @@ #include namespace times { + +/// Base class for correlation time distributions. +/// +/// Models the distribution of correlation times tau for individual walkers. +/// Each walker draws its own tau from the distribution during initialize(), +/// then waiting times between jumps are drawn from an exponential distribution +/// with mean tau_jump (Poisson process). +/// +/// Delta: all walkers have the same tau (no distribution). +/// LogNormal: tau is drawn from a log-normal distribution with width sigma. class BaseDistribution { public: virtual ~BaseDistribution() = default; @@ -21,19 +31,23 @@ public: [[nodiscard]] virtual std::unordered_map getParameters() const; + /// Draw a walker-specific correlation time from the distribution. virtual void initialize(std::mt19937_64 &rng) = 0; - [[nodiscard]] double tau_wait(std::mt19937_64 &rng) const; - [[nodiscard]] virtual std::unique_ptr clone() const = 0; + /// Draw a waiting time from the exponential distribution with mean m_tau_jump. + [[nodiscard]] double tau_wait(std::mt19937_64 &rng) const; + + [[nodiscard]] virtual std::unique_ptr clone() const = 0; [[nodiscard]] virtual std::string toString() const = 0; + /// Factory: create a distribution by registered name. static std::unique_ptr createFromInput(const std::string &input); protected: std::string m_name{"BaseDistribution"}; - double m_tau{1.}; - double m_tau_jump{1.}; + double m_tau{1.}; ///< Mean correlation time (center of distribution) + double m_tau_jump{1.}; ///< Walker-specific correlation time (drawn in initialize) }; } // namespace times diff --git a/src/utils/functions.cpp b/src/utils/functions.cpp index ae62913..110bafa 100644 --- a/src/utils/functions.cpp +++ b/src/utils/functions.cpp @@ -3,6 +3,9 @@ #include #include +/// Find index i such that x_ref[i] <= x < x_ref[i+1]. +/// When start > 0, scans linearly from the hint (efficient for sequential lookups). +/// When start == 0, uses binary search on large vectors. int nearest_index(const std::vector &x_ref, const double x, int start = 0) { const int last = static_cast(x_ref.size()) - 2; @@ -20,9 +23,7 @@ int nearest_index(const std::vector &x_ref, const double x, double lerp(const std::vector &x_ref, const std::vector &y_ref, const double x, const int i) { - /* - * Linear interpolation between two - */ + // Linear interpolation between points (x_ref[i], y_ref[i]) and (x_ref[i+1], y_ref[i+1]). const double x_left = x_ref[i]; const double y_left = y_ref[i]; const double x_right = x_ref[i + 1]; diff --git a/src/utils/registry.h b/src/utils/registry.h index 9127976..6316cc7 100644 --- a/src/utils/registry.h +++ b/src/utils/registry.h @@ -8,6 +8,14 @@ #include #include +/// Self-registering factory for polymorphic types. +/// +/// Models register themselves via static AutoRegister objects in their .cpp files: +/// static AutoRegister reg("RandomJump"); +/// +/// Important: libraries containing AutoRegister statics must be OBJECT libraries +/// (not STATIC) in CMake, otherwise the linker strips unreferenced translation +/// units and the registrations never execute. template class Registry { public: using Creator = std::function()>; @@ -21,6 +29,8 @@ public: entries()[name] = std::move(creator); } + /// Create an instance by name. Throws invalid_argument with available names + /// if the name is not found. static std::unique_ptr create(const std::string &name) { auto &map = entries(); auto it = map.find(name); @@ -39,6 +49,7 @@ public: } }; +/// Place a static instance in each model's .cpp to register it with the factory. template struct AutoRegister { explicit AutoRegister(const std::string &name) { Registry::add(name, []() { return std::make_unique(); });