From ef3cd31be615e87c157df8e142f2bca870e5a614 Mon Sep 17 00:00:00 2001 From: Dominik Demuth Date: Sun, 8 Mar 2026 11:28:02 +0100 Subject: [PATCH] Use unique_ptr instead of raw pointer --- src/main.cpp | 8 ++-- src/motions/base.cpp | 34 ++++++++-------- src/motions/base.h | 17 ++++---- src/motions/bimodalangle.cpp | 20 ++++++---- src/motions/bimodalangle.h | 9 +++-- src/motions/conewobble.cpp | 23 ++++++----- src/motions/conewobble.h | 14 +++---- src/motions/foursitejump.cpp | 26 ++++++------ src/motions/foursitejump.h | 10 ++--- src/motions/isosmallangle.cpp | 20 +++++----- src/motions/isosmallangle.h | 9 +++-- src/motions/nsiteconejump.cpp | 27 ++++++------- src/motions/nsiteconejump.h | 10 ++--- src/motions/random.cpp | 16 +++++--- src/motions/random.h | 10 ++--- src/motions/rjoac.cpp | 18 +++++---- src/motions/rjoac.h | 14 +++---- src/motions/sixsitejump.cpp | 28 ++++++------- src/motions/sixsitejump.h | 12 +++--- src/sims.cpp | 75 ++++++++++++++++------------------- src/sims.h | 39 +++++------------- src/times/base.cpp | 16 ++++---- src/times/base.h | 14 +++---- src/times/delta.cpp | 10 +++-- src/times/delta.h | 8 ++-- src/times/lognormal.cpp | 12 +++--- src/times/lognormal.h | 8 ++-- 27 files changed, 248 insertions(+), 259 deletions(-) diff --git a/src/main.cpp b/src/main.cpp index a249a8c..e415fb0 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -41,13 +41,13 @@ int main (const int argc, char *argv[]) { rng.seed(rd()); } - motions::BaseMotion *motion = motions::BaseMotion::createFromInput(args.motion_type, rng); - times::BaseDistribution *dist = times::BaseDistribution::createFromInput(args.distribution_type, rng); + auto motion = motions::BaseMotion::createFromInput(args.motion_type); + auto dist = times::BaseDistribution::createFromInput(args.distribution_type); if (args.spectrum) { - run_spectrum(parameter, args.optional, *motion, *dist); + run_spectrum(parameter, args.optional, *motion, *dist, rng); } if (args.ste) { - run_ste(parameter, args.optional, *motion, *dist); + run_ste(parameter, args.optional, *motion, *dist, rng); } return 0; diff --git a/src/motions/base.cpp b/src/motions/base.cpp index 96ee8f2..3de44df 100644 --- a/src/motions/base.cpp +++ b/src/motions/base.cpp @@ -11,19 +11,16 @@ #include "nsiteconejump.h" #include "sixsitejump.h" #include "rjoac.h" +#include "conewobble.h" #include namespace motions { - BaseMotion::BaseMotion(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.); - } + BaseMotion::BaseMotion(std::string name, const double delta, const double eta) : m_name(std::move(name)), m_delta(delta), m_eta(eta) {} - BaseMotion::BaseMotion(std::string name, std::mt19937_64& rng) : m_name(std::move(name)), m_rng(rng) { - m_uni_dist = std::uniform_real_distribution(0., 1.); - } + BaseMotion::BaseMotion(std::string name) : m_name(std::move(name)) {} double BaseMotion::omega_q(const double cos_theta, const double phi) const { const double cos_theta_square = cos_theta * cos_theta; @@ -38,34 +35,37 @@ namespace motions { return omega_q(cos_theta, phi); } - coordinates::SphericalPos BaseMotion::draw_position() { - const double cos_theta = 1 - 2 * m_uni_dist(m_rng); - const double phi = 2.0 * M_PI * m_uni_dist(m_rng); + coordinates::SphericalPos BaseMotion::draw_position(std::mt19937_64& rng) { + const double cos_theta = 1 - 2 * m_uni_dist(rng); + const double phi = 2.0 * M_PI * m_uni_dist(rng); return {cos_theta, phi}; } - BaseMotion* BaseMotion::createFromInput(const std::string& input, std::mt19937_64& rng) { + std::unique_ptr BaseMotion::createFromInput(const std::string& input) { if (input == "FourSiteTetrahedral") - return new FourSiteTetrahedron(rng); + return std::make_unique(); if (input == "SixSiteOctahedralC3") - return new SixSiteOctahedronC3(rng); + return std::make_unique(); if (input == "IsotropicAngle") - return new SmallAngle(rng); + return std::make_unique(); if (input == "RandomJump") - return new RandomJump(rng); + return std::make_unique(); if (input == "BimodalAngle") - return new BimodalAngle(rng); + return std::make_unique(); if (input == "NSiteConeJump") - return new NSiteJumpOnCone(rng); + return std::make_unique(); if (input == "RandomJumpOnCone") - return new RandomJumpOnCone(rng); + return std::make_unique(); + + if (input == "ConeWobble") + return std::make_unique(); throw std::invalid_argument("Invalid input " + input); } diff --git a/src/motions/base.h b/src/motions/base.h index 68f6961..cc37f23 100644 --- a/src/motions/base.h +++ b/src/motions/base.h @@ -3,6 +3,7 @@ #include "coordinates.h" +#include #include #include @@ -11,15 +12,16 @@ namespace motions { public: virtual ~BaseMotion() = default; - BaseMotion(std::string, double, double, std::mt19937_64&); - explicit BaseMotion(std::string, std::mt19937_64&); + BaseMotion(std::string, double, double); + explicit BaseMotion(std::string); - coordinates::SphericalPos draw_position(); + coordinates::SphericalPos draw_position(std::mt19937_64& rng); [[nodiscard]] double omega_q(double, double) const; [[nodiscard]] double omega_q(const coordinates::SphericalPos&) const; - virtual void initialize() = 0; - virtual double jump() = 0; + virtual void initialize(std::mt19937_64& rng) = 0; + virtual double jump(std::mt19937_64& rng) = 0; + [[nodiscard]] virtual std::unique_ptr clone() const = 0; virtual void setParameters(const std::unordered_map&); [[nodiscard]] virtual std::unordered_map getParameters() const; @@ -32,14 +34,13 @@ namespace motions { [[nodiscard]] virtual std::string toString() const = 0; - static BaseMotion* createFromInput(const std::string& input, std::mt19937_64& rng); + static std::unique_ptr createFromInput(const std::string& input); protected: std::string m_name{"BaseMotion"}; double m_delta{1.}; double m_eta{0.}; - std::mt19937_64& m_rng; - std::uniform_real_distribution<> m_uni_dist; + std::uniform_real_distribution<> m_uni_dist{0., 1.}; double m_initial_omega{0.}; }; diff --git a/src/motions/bimodalangle.cpp b/src/motions/bimodalangle.cpp index be64719..ce33de5 100644 --- a/src/motions/bimodalangle.cpp +++ b/src/motions/bimodalangle.cpp @@ -4,26 +4,30 @@ #include "coordinates.h" namespace motions { - BimodalAngle::BimodalAngle(const double delta, const double eta, const double angle1, const double angle2, const double prob, std::mt19937_64 &rng) : - BaseMotion(std::string("BimodalAngle"), delta, eta, rng), + BimodalAngle::BimodalAngle(const double delta, const double eta, const double angle1, const double angle2, const double prob) : + BaseMotion(std::string("BimodalAngle"), delta, eta), m_angle1(angle1 * M_PI / 180.0), m_angle2(angle2 * M_PI / 180.0), m_prob(prob) {} - BimodalAngle::BimodalAngle(std::mt19937_64 &rng) : BaseMotion(std::string("BimodalAngle"), rng) {} + BimodalAngle::BimodalAngle() : BaseMotion(std::string("BimodalAngle")) {} - void BimodalAngle::initialize() { - m_prev_pos = draw_position(); + void BimodalAngle::initialize(std::mt19937_64& rng) { + m_prev_pos = draw_position(rng); m_initial_omega = omega_q(m_prev_pos); } - double BimodalAngle::jump() { - const double angle = m_uni_dist(m_rng) < m_prob ? m_angle1 : m_angle2; - const double gamma{2 * M_PI * m_uni_dist(m_rng)}; + double BimodalAngle::jump(std::mt19937_64& rng) { + const double angle = m_uni_dist(rng) < m_prob ? m_angle1 : m_angle2; + const double gamma{2 * M_PI * m_uni_dist(rng)}; m_prev_pos = rotate(m_prev_pos, angle, gamma); return omega_q(m_prev_pos); } + std::unique_ptr BimodalAngle::clone() const { + return std::make_unique(*this); + } + void BimodalAngle::setParameters(const std::unordered_map ¶meter) { BaseMotion::setParameters(parameter); diff --git a/src/motions/bimodalangle.h b/src/motions/bimodalangle.h index 6b05ace..491fb76 100644 --- a/src/motions/bimodalangle.h +++ b/src/motions/bimodalangle.h @@ -8,11 +8,12 @@ namespace motions { class BimodalAngle final : public BaseMotion { public: - BimodalAngle(double, double, double, double, double, std::mt19937_64& ); - explicit BimodalAngle(std::mt19937_64&); + BimodalAngle(double, double, double, double, double); + BimodalAngle(); - void initialize() override; - double jump() override; + void initialize(std::mt19937_64& rng) override; + double jump(std::mt19937_64& rng) override; + [[nodiscard]] std::unique_ptr clone() const override; void setParameters(const std::unordered_map &) override; [[nodiscard]] std::unordered_map getParameters() const override; [[nodiscard]] std::string toString() const override; diff --git a/src/motions/conewobble.cpp b/src/motions/conewobble.cpp index 660d7c1..db0af4c 100644 --- a/src/motions/conewobble.cpp +++ b/src/motions/conewobble.cpp @@ -2,23 +2,24 @@ #include "conewobble.h" #include "coordinates.h" -#include -#include - namespace motions { - WobbleCone::WobbleCone(const double delta, const double eta, const double chi, std::mt19937_64 &rng) : BaseMotion("Wobble in Cone", delta, eta, rng), m_angle(chi) {} - WobbleCone::WobbleCone(std::mt19937_64 &rng) : BaseMotion("Wobble in Cone", rng) {} + WobbleCone::WobbleCone(const double delta, const double eta, const double chi) : BaseMotion("Wobble in Cone", delta, eta), m_angle(chi) {} + WobbleCone::WobbleCone() : BaseMotion("Wobble in Cone") {} - void WobbleCone::initialize() { - m_axis = draw_position(); + void WobbleCone::initialize(std::mt19937_64& rng) { + m_axis = draw_position(rng); } - double WobbleCone::jump() { - const double real_angle = m_uni_dist(m_rng) * m_angle; - const double phi = 2 * M_PI * m_uni_dist(m_rng); + double WobbleCone::jump(std::mt19937_64& rng) { + const double real_angle = m_uni_dist(rng) * m_angle; + const double phi = 2 * M_PI * m_uni_dist(rng); return omega_q(rotate(m_axis, real_angle, phi)); } + std::unique_ptr WobbleCone::clone() const { + return std::make_unique(*this); + } + void WobbleCone::setParameters(const std::unordered_map ¶meters) { BaseMotion::setParameters(parameters); m_angle = parameters.at("angle"); @@ -34,4 +35,4 @@ namespace motions { std::string WobbleCone::toString() const { return std::string("ConeWobble/angle=") + std::to_string(m_angle); } -} \ No newline at end of file +} diff --git a/src/motions/conewobble.h b/src/motions/conewobble.h index 715be43..1c42a78 100644 --- a/src/motions/conewobble.h +++ b/src/motions/conewobble.h @@ -4,22 +4,18 @@ #include "base.h" #include "coordinates.h" -#include - namespace motions { class WobbleCone final: public BaseMotion { public: - WobbleCone(double, double, double, std::mt19937_64&); - explicit WobbleCone(std::mt19937_64&); + WobbleCone(double, double, double); + WobbleCone(); - void initialize() override; - - double jump() override; + void initialize(std::mt19937_64& rng) override; + double jump(std::mt19937_64& rng) override; + [[nodiscard]] std::unique_ptr clone() const override; void setParameters(const std::unordered_map &) override; - [[nodiscard]] std::unordered_map getParameters() const override; - [[nodiscard]] std::string toString() const override; private: diff --git a/src/motions/foursitejump.cpp b/src/motions/foursitejump.cpp index 40fc5cf..464b60b 100644 --- a/src/motions/foursitejump.cpp +++ b/src/motions/foursitejump.cpp @@ -1,35 +1,37 @@ #include "foursitejump.h" #include "coordinates.h" -#include - namespace motions { - FourSiteTetrahedron::FourSiteTetrahedron(const double delta, const double eta, std::mt19937_64& rng) : - BaseMotion(std::string{"FourSiteTetrahedral"}, delta, eta, rng) {} + FourSiteTetrahedron::FourSiteTetrahedron(const double delta, const double eta) : + BaseMotion(std::string{"FourSiteTetrahedral"}, delta, eta) {} - FourSiteTetrahedron::FourSiteTetrahedron(std::mt19937_64& rng) : BaseMotion(std::string{"FourSiteTetrahedral"}, rng) {} + FourSiteTetrahedron::FourSiteTetrahedron() : BaseMotion(std::string{"FourSiteTetrahedral"}) {} - void FourSiteTetrahedron::initialize() { - const auto pos = draw_position(); + void FourSiteTetrahedron::initialize(std::mt19937_64& rng) { + const auto pos = draw_position(rng); m_corners[0] = omega_q(pos); - const double alpha = 2. * M_PI * m_uni_dist(m_rng); + const double alpha = 2. * M_PI * m_uni_dist(rng); for (int i = 1; i<4; i++) { auto corner_pos = rotate(pos, m_beta, alpha + (i-1) * 2*M_PI/3.); m_corners[i] = omega_q(corner_pos); } - m_initial_omega = FourSiteTetrahedron::jump(); + m_initial_omega = FourSiteTetrahedron::jump(rng); } - double FourSiteTetrahedron::jump() { - m_corner_idx += m_chooser(m_rng); + double FourSiteTetrahedron::jump(std::mt19937_64& rng) { + m_corner_idx += m_chooser(rng); m_corner_idx %= 4; return m_corners[m_corner_idx]; } + std::unique_ptr FourSiteTetrahedron::clone() const { + return std::make_unique(*this); + } + std::string FourSiteTetrahedron::toString() const { return {"FourSiteTetrahedral"}; } -} \ No newline at end of file +} diff --git a/src/motions/foursitejump.h b/src/motions/foursitejump.h index 7ae8538..19bc0db 100644 --- a/src/motions/foursitejump.h +++ b/src/motions/foursitejump.h @@ -3,18 +3,18 @@ #include "base.h" -#include #include #include namespace motions { class FourSiteTetrahedron final : public BaseMotion { public: - FourSiteTetrahedron(double, double, std::mt19937_64&); - explicit FourSiteTetrahedron(std::mt19937_64&); + FourSiteTetrahedron(double, double); + FourSiteTetrahedron(); - void initialize() override; - double jump() override; + void initialize(std::mt19937_64& rng) override; + double jump(std::mt19937_64& rng) override; + [[nodiscard]] std::unique_ptr clone() const override; [[nodiscard]] std::string toString() const override; diff --git a/src/motions/isosmallangle.cpp b/src/motions/isosmallangle.cpp index be6ebeb..94beb43 100644 --- a/src/motions/isosmallangle.cpp +++ b/src/motions/isosmallangle.cpp @@ -1,26 +1,28 @@ #include "isosmallangle.h" #include "coordinates.h" -#include - namespace motions { - SmallAngle::SmallAngle(const double delta, const double eta, const double chi, std::mt19937_64 &rng) : - BaseMotion(std::string("IsotropicAngle"), delta, eta, rng), m_chi(chi * M_PI / 180.0) {} - SmallAngle::SmallAngle(std::mt19937_64 &rng) : BaseMotion(std::string("IsotropicAngle"), rng) {} + SmallAngle::SmallAngle(const double delta, const double eta, const double chi) : + BaseMotion(std::string("IsotropicAngle"), delta, eta), m_chi(chi * M_PI / 180.0) {} + SmallAngle::SmallAngle() : BaseMotion(std::string("IsotropicAngle")) {} - void SmallAngle::initialize() { - m_prev_pos = draw_position(); + void SmallAngle::initialize(std::mt19937_64& rng) { + m_prev_pos = draw_position(rng); m_initial_omega = omega_q(m_prev_pos); } - double SmallAngle::jump() { - const double gamma{2 * M_PI * m_uni_dist(m_rng)}; + double SmallAngle::jump(std::mt19937_64& rng) { + const double gamma{2 * M_PI * m_uni_dist(rng)}; m_prev_pos = rotate(m_prev_pos, m_chi, gamma); return omega_q(m_prev_pos); } + std::unique_ptr SmallAngle::clone() const { + return std::make_unique(*this); + } + void SmallAngle::setParameters(const std::unordered_map ¶meters) { m_chi = parameters.at("angle") * M_PI / 180.0; BaseMotion::setParameters(parameters); diff --git a/src/motions/isosmallangle.h b/src/motions/isosmallangle.h index 0e04f16..ed14e09 100644 --- a/src/motions/isosmallangle.h +++ b/src/motions/isosmallangle.h @@ -8,11 +8,12 @@ namespace motions { class SmallAngle final : public BaseMotion { public: - SmallAngle(double, double, double, std::mt19937_64& ); - explicit SmallAngle(std::mt19937_64&); + SmallAngle(double, double, double); + SmallAngle(); - void initialize() override; - double jump() override; + void initialize(std::mt19937_64& rng) override; + double jump(std::mt19937_64& rng) override; + [[nodiscard]] std::unique_ptr clone() const override; void setParameters(const std::unordered_map &) override; [[nodiscard]] std::unordered_map getParameters() const override; [[nodiscard]] std::string toString() const override; diff --git a/src/motions/nsiteconejump.cpp b/src/motions/nsiteconejump.cpp index 56976dd..b737d4f 100644 --- a/src/motions/nsiteconejump.cpp +++ b/src/motions/nsiteconejump.cpp @@ -1,43 +1,42 @@ -// -// Created by dominik on 12/14/24. -// #include "nsiteconejump.h" #include "coordinates.h" #include -#include -#include #include namespace motions { - NSiteJumpOnCone::NSiteJumpOnCone(const double delta, const double eta, const int num_sites, const double chi, std::mt19937_64 &rng) : - BaseMotion("NSiteJumpOnCone", delta, eta, rng), + NSiteJumpOnCone::NSiteJumpOnCone(const double delta, const double eta, const int num_sites, const double chi) : + BaseMotion("NSiteJumpOnCone", delta, eta), m_num_sites(num_sites), m_chi(chi*M_PI/180.) {} - NSiteJumpOnCone::NSiteJumpOnCone(std::mt19937_64 &rng) : BaseMotion("NSiteJumpOnCone", rng) { } + NSiteJumpOnCone::NSiteJumpOnCone() : BaseMotion("NSiteJumpOnCone") { } - void NSiteJumpOnCone::initialize() { + void NSiteJumpOnCone::initialize(std::mt19937_64& rng) { m_sites = std::vector(m_num_sites); m_chooser = std::uniform_int_distribution<>{1, m_num_sites - 1}; - m_axis = draw_position(); + m_axis = draw_position(rng); - const double alpha = m_uni_dist(m_rng) * 2 * M_PI; + const double alpha = m_uni_dist(rng) * 2 * M_PI; const double steps = 2*M_PI / m_num_sites; for (int i = 0; i < m_num_sites; i++) { m_sites[i] = omega_q(rotate(m_axis, m_chi, i * steps + alpha)); } } - double NSiteJumpOnCone::jump() { - m_cone_idx += m_chooser(m_rng); + double NSiteJumpOnCone::jump(std::mt19937_64& rng) { + m_cone_idx += m_chooser(rng); m_cone_idx %= m_num_sites; return m_sites[m_cone_idx]; } + std::unique_ptr NSiteJumpOnCone::clone() const { + return std::make_unique(*this); + } + void NSiteJumpOnCone::setParameters(const std::unordered_map ¶meters) { BaseMotion::setParameters(parameters); m_num_sites = static_cast(parameters.at("num_sites")); @@ -56,4 +55,4 @@ namespace motions { return std::to_string(m_num_sites) + "SiteJumpOnCone/angle=" + std::to_string(m_chi*180/M_PI); } -} // motions \ No newline at end of file +} // motions diff --git a/src/motions/nsiteconejump.h b/src/motions/nsiteconejump.h index 116f620..3d0f3d7 100644 --- a/src/motions/nsiteconejump.h +++ b/src/motions/nsiteconejump.h @@ -3,17 +3,17 @@ #include "base.h" -#include #include namespace motions { class NSiteJumpOnCone final : public BaseMotion { public: - NSiteJumpOnCone(double, double, int, double, std::mt19937_64&); - explicit NSiteJumpOnCone(std::mt19937_64&); + NSiteJumpOnCone(double, double, int, double); + NSiteJumpOnCone(); - void initialize() override; - double jump() override; + void initialize(std::mt19937_64& rng) override; + double jump(std::mt19937_64& rng) override; + [[nodiscard]] std::unique_ptr clone() const override; [[nodiscard]] std::string toString() const override; diff --git a/src/motions/random.cpp b/src/motions/random.cpp index eb34481..f4b786f 100644 --- a/src/motions/random.cpp +++ b/src/motions/random.cpp @@ -2,19 +2,23 @@ #include "random.h" namespace motions { - RandomJump::RandomJump(const double delta, const double eta, std::mt19937_64 &rng) : BaseMotion(std::string("RandomJump"), delta, eta, rng) {} + RandomJump::RandomJump(const double delta, const double eta) : BaseMotion(std::string("RandomJump"), delta, eta) {} - RandomJump::RandomJump(std::mt19937_64 &rng) : BaseMotion(std::string("RandomJump"), rng) {} + RandomJump::RandomJump() : BaseMotion(std::string("RandomJump")) {} std::string RandomJump::toString() const { return {"RandomJump"}; } - void RandomJump::initialize() { - m_initial_omega = RandomJump::jump(); + std::unique_ptr RandomJump::clone() const { + return std::make_unique(*this); } - double RandomJump::jump() { - return omega_q(draw_position()); + void RandomJump::initialize(std::mt19937_64& rng) { + m_initial_omega = RandomJump::jump(rng); + } + + double RandomJump::jump(std::mt19937_64& rng) { + return omega_q(draw_position(rng)); } } diff --git a/src/motions/random.h b/src/motions/random.h index 3eb465b..ea86246 100644 --- a/src/motions/random.h +++ b/src/motions/random.h @@ -3,18 +3,18 @@ #include "base.h" -#include namespace motions { class RandomJump final : public BaseMotion { public: - RandomJump(double, double, std::mt19937_64&); - explicit RandomJump(std::mt19937_64&); + RandomJump(double, double); + RandomJump(); [[nodiscard]] std::string toString() const override; + [[nodiscard]] std::unique_ptr clone() const override; - void initialize() override; - double jump() override; + void initialize(std::mt19937_64& rng) override; + double jump(std::mt19937_64& rng) override; }; } diff --git a/src/motions/rjoac.cpp b/src/motions/rjoac.cpp index 098e848..57402fd 100644 --- a/src/motions/rjoac.cpp +++ b/src/motions/rjoac.cpp @@ -2,18 +2,22 @@ #include "rjoac.h" namespace motions { - RandomJumpOnCone::RandomJumpOnCone(const double delta, const double eta, const double chi, std::mt19937_64 &rng) : BaseMotion("RJ on a Cone", delta, eta, rng), m_angle(chi) {} - RandomJumpOnCone::RandomJumpOnCone(std::mt19937_64 &rng) : BaseMotion("RJ on a Cone", rng) {} + RandomJumpOnCone::RandomJumpOnCone(const double delta, const double eta, const double chi) : BaseMotion("RJ on a Cone", delta, eta), m_angle(chi) {} + RandomJumpOnCone::RandomJumpOnCone() : BaseMotion("RJ on a Cone") {} - void RandomJumpOnCone::initialize() { - m_axis = draw_position(); + void RandomJumpOnCone::initialize(std::mt19937_64& rng) { + m_axis = draw_position(rng); } - double RandomJumpOnCone::jump() { - const double phi = 2 * M_PI * m_uni_dist(m_rng); + double RandomJumpOnCone::jump(std::mt19937_64& rng) { + const double phi = 2 * M_PI * m_uni_dist(rng); return omega_q(rotate(m_axis, m_angle, phi)); } + std::unique_ptr RandomJumpOnCone::clone() const { + return std::make_unique(*this); + } + void RandomJumpOnCone::setParameters(const std::unordered_map ¶meters) { BaseMotion::setParameters(parameters); m_angle = parameters.at("angle"); @@ -29,4 +33,4 @@ namespace motions { std::string RandomJumpOnCone::toString() const { return std::string("RandomJumpOnCone/angle=") + std::to_string(m_angle); } -} \ No newline at end of file +} diff --git a/src/motions/rjoac.h b/src/motions/rjoac.h index d96a7ac..19392e5 100644 --- a/src/motions/rjoac.h +++ b/src/motions/rjoac.h @@ -4,22 +4,18 @@ #include "base.h" #include "coordinates.h" -#include - namespace motions { class RandomJumpOnCone final: public BaseMotion { public: - RandomJumpOnCone(double, double, double, std::mt19937_64&); - explicit RandomJumpOnCone(std::mt19937_64&); + RandomJumpOnCone(double, double, double); + RandomJumpOnCone(); - void initialize() override; - - double jump() override; + void initialize(std::mt19937_64& rng) override; + double jump(std::mt19937_64& rng) override; + [[nodiscard]] std::unique_ptr clone() const override; void setParameters(const std::unordered_map &) override; - [[nodiscard]] std::unordered_map getParameters() const override; - [[nodiscard]] std::string toString() const override; private: diff --git a/src/motions/sixsitejump.cpp b/src/motions/sixsitejump.cpp index c34bf7d..0f859ae 100644 --- a/src/motions/sixsitejump.cpp +++ b/src/motions/sixsitejump.cpp @@ -1,21 +1,17 @@ #include "sixsitejump.h" -#include -#include - #include "coordinates.h" namespace motions { - SixSiteOctahedronC3::SixSiteOctahedronC3(const double delta, const double eta, const double chi, std::mt19937_64& rng) : - m_chi{chi*M_PI/180.}, - BaseMotion(std::string{"SixSiteOctahedral"}, delta, eta, rng) {} + SixSiteOctahedronC3::SixSiteOctahedronC3(const double delta, const double eta, const double chi) : + BaseMotion(std::string{"SixSiteOctahedral"}, delta, eta), + m_chi{chi*M_PI/180.} {} - SixSiteOctahedronC3::SixSiteOctahedronC3(std::mt19937_64& rng) : BaseMotion(std::string{"SixSiteOctahedralC3"}, rng) {} + SixSiteOctahedronC3::SixSiteOctahedronC3() : BaseMotion(std::string{"SixSiteOctahedralC3"}) {} - void SixSiteOctahedronC3::initialize() { - const coordinates::SphericalPos c3_axis = draw_position(); - const auto [x, y, z] = spherical_to_xyz(c3_axis); - const double alpha = 2. * M_PI * m_uni_dist(m_rng); + void SixSiteOctahedronC3::initialize(std::mt19937_64& rng) { + const coordinates::SphericalPos c3_axis = draw_position(rng); + const double alpha = 2. * M_PI * m_uni_dist(rng); const double m_chi_opposite = M_PI - m_chi; @@ -24,17 +20,21 @@ namespace motions { m_corners[2*i+1] = omega_q(rotate(c3_axis, m_chi_opposite, alpha + i * 2./3.*M_PI + M_PI/3.)); } - m_initial_omega = SixSiteOctahedronC3::jump(); + m_initial_omega = SixSiteOctahedronC3::jump(rng); } - double SixSiteOctahedronC3::jump() { - m_corner_idx += m_chooser(m_rng); + double SixSiteOctahedronC3::jump(std::mt19937_64& rng) { + m_corner_idx += m_chooser(rng); m_corner_idx %= 6; return m_corners[m_corner_idx]; } + std::unique_ptr SixSiteOctahedronC3::clone() const { + return std::make_unique(*this); + } + std::string SixSiteOctahedronC3::toString() const { return {"SixSiteOctahedral/angle=" + std::to_string(m_chi / M_PI * 180.)}; } diff --git a/src/motions/sixsitejump.h b/src/motions/sixsitejump.h index 8e203f6..314de74 100644 --- a/src/motions/sixsitejump.h +++ b/src/motions/sixsitejump.h @@ -4,23 +4,23 @@ #include "base.h" -#include #include #include namespace motions { class SixSiteOctahedronC3 final : public BaseMotion { public: - SixSiteOctahedronC3(double, double, double, std::mt19937_64&); - explicit SixSiteOctahedronC3(std::mt19937_64&); + SixSiteOctahedronC3(double, double, double); + SixSiteOctahedronC3(); - void initialize() override; - double jump() override; + void initialize(std::mt19937_64& rng) override; + double jump(std::mt19937_64& rng) override; + [[nodiscard]] std::unique_ptr clone() const override; [[nodiscard]] std::string toString() const override; private: - const double m_chi{0.95531661812450927816385710251575775424341469501000549095969812932191204590}; // 54.74 deg + double m_chi{0.95531661812450927816385710251575775424341469501000549095969812932191204590}; // 54.74 deg std::array m_corners{}; int m_corner_idx{0}; diff --git a/src/sims.cpp b/src/sims.cpp index 3f14a59..f07e824 100644 --- a/src/sims.cpp +++ b/src/sims.cpp @@ -19,7 +19,8 @@ void run_spectrum( std::unordered_map& parameter, std::unordered_map optional, motions::BaseMotion& motion, - times::BaseDistribution& dist + times::BaseDistribution& dist, + std::mt19937_64& rng ) { const int num_walker = static_cast(parameter["num_walker"]); @@ -47,23 +48,19 @@ void run_spectrum( // let the walker walk 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, traj_omega); + auto traj = make_trajectory(motion, dist, tmax, rng); for (auto& [t_echo_j, fid_j] : fid_dict) { // get phase at echo pulse - int current_pos = nearest_index(traj_time, t_echo_j, 0); - const double phase_techo = lerp(traj_time, traj_phase, t_echo_j, current_pos); + int current_pos = nearest_index(traj.time, t_echo_j, 0); + const double phase_techo = lerp(traj.time, traj.phase, t_echo_j, current_pos); for (int acq_idx = 0; acq_idx < num_acq; acq_idx++) { const double real_time = t_fid[acq_idx] + 2 * t_echo_j; - current_pos = nearest_index(traj_time, real_time, current_pos); - const double phase_acq = lerp(traj_time, traj_phase, real_time, current_pos); + current_pos = nearest_index(traj.time, real_time, current_pos); + const double phase_acq = lerp(traj.time, traj.phase, real_time, current_pos); fid_j[acq_idx] += std::cos(phase_acq - 2 * phase_techo) / num_walker; } @@ -84,7 +81,8 @@ void run_ste( std::unordered_map& parameter, std::unordered_map optional, motions::BaseMotion& motion, - times::BaseDistribution& dist + times::BaseDistribution& dist, + std::mt19937_64& rng ) { const int num_walker = static_cast(parameter[std::string("num_walker")]); @@ -116,17 +114,13 @@ void run_ste( // let the walker walk 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, traj_omega); + auto traj = make_trajectory(motion, dist, tmax, rng); 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; + f2_pos = nearest_index(traj.time, t_mix_f2, f2_pos); + f2[f2_idx] += traj.omega[f2_pos] * motion.getInitOmega() / num_walker; } @@ -135,26 +129,26 @@ void run_ste( auto& ss_j = ss_dict[t_evo_j]; // get phase at beginning of mixing time - int current_pos = nearest_index(traj_time, t_evo_j, 0); - const double dephased = lerp(traj_time, traj_phase, t_evo_j, current_pos); + int current_pos = nearest_index(traj.time, t_evo_j, 0); + const double dephased = lerp(traj.time, traj.phase, t_evo_j, current_pos); const double cc_tevo = std::cos(dephased); const double ss_tevo = std::sin(dephased); for (int mix_idx = 0; mix_idx < num_mix_times; mix_idx++) { // get phase at end of mixing time const double time_end_mix = mixing_times[mix_idx] + t_evo_j; - current_pos = nearest_index(traj_time, time_end_mix, current_pos); - const double phase_mix_end = lerp(traj_time, traj_phase, time_end_mix, current_pos); + current_pos = nearest_index(traj.time, time_end_mix, current_pos); + const double phase_mix_end = lerp(traj.time, traj.phase, time_end_mix, current_pos); // get phase at position of 4th pulse const double time_pulse4 = time_end_mix + tpulse4; - current_pos = nearest_index(traj_time, time_pulse4, current_pos); - const double phase_4pulse = lerp(traj_time, traj_phase, time_pulse4, current_pos); + current_pos = nearest_index(traj.time, time_pulse4, current_pos); + const double phase_4pulse = lerp(traj.time, traj.phase, time_pulse4, current_pos); // get phase at echo position const double time_echo = time_pulse4 + tpulse4 + t_evo_j; - current_pos = nearest_index(traj_time, time_echo, current_pos); - double rephased = lerp(traj_time, traj_phase, time_echo, current_pos) + phase_mix_end - 2*phase_4pulse; + current_pos = nearest_index(traj.time, time_echo, current_pos); + double rephased = lerp(traj.time, traj.phase, time_echo, current_pos) + phase_mix_end - 2*phase_4pulse; cc_j[mix_idx] += cc_tevo * std::cos(rephased) / num_walker; ss_j[mix_idx] += ss_tevo * std::sin(rephased) / num_walker; @@ -174,37 +168,38 @@ void run_ste( } -void make_trajectory( +Trajectory make_trajectory( motions::BaseMotion& motion, times::BaseDistribution& dist, const double t_max, - std::vector& out_time, - std::vector& out_phase, - std::vector& out_omega + std::mt19937_64& rng ) { // Starting position double t_passed = 0; double phase = 0; - motion.initialize(); - dist.initialize(); + motion.initialize(rng); + dist.initialize(rng); double omega = motion.getInitOmega(); - out_time.emplace_back(t_passed); - out_phase.emplace_back(phase); - out_omega.emplace_back(omega); + Trajectory traj; + traj.time.emplace_back(t_passed); + traj.phase.emplace_back(phase); + traj.omega.emplace_back(omega); while (t_passed < t_max) { - const double t = dist.tau_wait(); + const double t = dist.tau_wait(rng); t_passed += t; - omega = motion.jump(); + omega = motion.jump(rng); phase += omega * t; - out_time.emplace_back(t_passed); - out_phase.emplace_back(phase); - out_omega.emplace_back(omega); + traj.time.emplace_back(t_passed); + traj.phase.emplace_back(phase); + traj.omega.emplace_back(omega); } + + return traj; } diff --git a/src/sims.h b/src/sims.h index 8414185..1c5fb4b 100644 --- a/src/sims.h +++ b/src/sims.h @@ -7,38 +7,19 @@ #include #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, std::unordered_map optional, motions::BaseMotion& motion, times::BaseDistribution& dist); +struct Trajectory { + std::vector time; + std::vector phase; + std::vector omega; +}; -/** - * @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, std::unordered_map optional, motions::BaseMotion& motion, times::BaseDistribution& dist); +void run_spectrum(std::unordered_map& parameter, std::unordered_map optional, motions::BaseMotion& motion, times::BaseDistribution& dist, std::mt19937_64& rng); -/** - * @brief Create trajectory of a single walker - * - * @param motion Motion model - * @param dist Distribution of correlation times - * @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(motions::BaseMotion& motion, times::BaseDistribution& dist, double t_max, std::vector& out_time, std::vector& out_phase, std::vector& out_omega); +void run_ste(std::unordered_map& parameter, std::unordered_map optional, motions::BaseMotion& motion, times::BaseDistribution& dist, std::mt19937_64& rng); + +Trajectory make_trajectory(motions::BaseMotion& motion, times::BaseDistribution& dist, 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/times/base.cpp b/src/times/base.cpp index 500e651..f93e614 100644 --- a/src/times/base.cpp +++ b/src/times/base.cpp @@ -5,12 +5,12 @@ #include namespace times { - BaseDistribution::BaseDistribution(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) {} + BaseDistribution::BaseDistribution(std::string name, const double tau) : m_name(std::move(name)), m_tau(tau), m_tau_jump(tau) {} - BaseDistribution::BaseDistribution(std::string name, std::mt19937_64 &rng) : m_name(std::move(name)), m_rng(rng) {} + BaseDistribution::BaseDistribution(std::string name) : m_name(std::move(name)) {} - double BaseDistribution::tau_wait() const { - return std::exponential_distribution(1./m_tau_jump)(m_rng); + double BaseDistribution::tau_wait(std::mt19937_64& rng) const { + return std::exponential_distribution(1./m_tau_jump)(rng); } void BaseDistribution::setParameters(const std::unordered_map ¶meters) { @@ -24,13 +24,13 @@ namespace times { } - BaseDistribution* BaseDistribution::createFromInput(const std::string& input, std::mt19937_64& rng) { + std::unique_ptr BaseDistribution::createFromInput(const std::string& input) { if (input == "Delta") - return new DeltaDistribution(rng); + return std::make_unique(); if (input == "LogNormal") - return new LogNormalDistribution(rng); + return std::make_unique(); 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 index f909e5f..1f23c68 100644 --- a/src/times/base.h +++ b/src/times/base.h @@ -1,6 +1,7 @@ #ifndef RWSIM_TIMESBASE_H #define RWSIM_TIMESBASE_H +#include #include #include @@ -9,8 +10,8 @@ namespace times { public: virtual ~BaseDistribution() = default; - BaseDistribution(std::string, double, std::mt19937_64&); - explicit BaseDistribution(std::string, std::mt19937_64&); + BaseDistribution(std::string, double); + explicit BaseDistribution(std::string); [[nodiscard]] double getTau() const { return m_tau; } void setTau(const double tau) { m_tau = tau; } @@ -19,19 +20,18 @@ namespace times { virtual void setParameters(const std::unordered_map&); [[nodiscard]] virtual std::unordered_map getParameters() const; - virtual void initialize() = 0; - virtual void draw_tau() = 0; - [[nodiscard]] double tau_wait() const; + 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; [[nodiscard]] virtual std::string toString() const = 0; - static BaseDistribution* createFromInput(const std::string& input, std::mt19937_64& rng); + static std::unique_ptr createFromInput(const std::string& input); protected: std::string m_name{"BaseDistribution"}; double m_tau{1.}; double m_tau_jump{1.}; - std::mt19937_64& m_rng; }; } diff --git a/src/times/delta.cpp b/src/times/delta.cpp index 8f59508..e2e8a8f 100644 --- a/src/times/delta.cpp +++ b/src/times/delta.cpp @@ -1,14 +1,16 @@ #include "delta.h" namespace times { - DeltaDistribution::DeltaDistribution(const double tau, std::mt19937_64& rng) : BaseDistribution(std::string("Delta"), tau, rng) {} - DeltaDistribution::DeltaDistribution(std::mt19937_64& rng) : BaseDistribution(std::string("Delta"), rng) {} + DeltaDistribution::DeltaDistribution(const double tau) : BaseDistribution(std::string("Delta"), tau) {} + DeltaDistribution::DeltaDistribution() : BaseDistribution(std::string("Delta")) {} - void DeltaDistribution::initialize() { + void DeltaDistribution::initialize(std::mt19937_64&) { m_tau_jump = m_tau; } - void DeltaDistribution::draw_tau() {} + std::unique_ptr DeltaDistribution::clone() const { + return std::make_unique(*this); + } std::string DeltaDistribution::toString() const { return {"Delta/tau=" + std::to_string(m_tau)}; diff --git a/src/times/delta.h b/src/times/delta.h index 33b527f..30a57b5 100644 --- a/src/times/delta.h +++ b/src/times/delta.h @@ -6,11 +6,11 @@ namespace times { class DeltaDistribution final : public BaseDistribution { public: - DeltaDistribution(double, std::mt19937_64&); - explicit DeltaDistribution(std::mt19937_64 &rng); + explicit DeltaDistribution(double); + DeltaDistribution(); - void initialize() override; - void draw_tau() override; + void initialize(std::mt19937_64& rng) override; + [[nodiscard]] std::unique_ptr clone() const override; [[nodiscard]] std::string toString() const override; }; diff --git a/src/times/lognormal.cpp b/src/times/lognormal.cpp index 023e1ea..1c15fd8 100644 --- a/src/times/lognormal.cpp +++ b/src/times/lognormal.cpp @@ -2,8 +2,8 @@ #include namespace times { - LogNormalDistribution::LogNormalDistribution(const double tau, const double sigma, std::mt19937_64& rng) : BaseDistribution(std::string("LogNormal"), tau, rng), m_sigma(sigma), m_distribution(std::log(tau), sigma) {} - LogNormalDistribution::LogNormalDistribution(std::mt19937_64& rng) : BaseDistribution(std::string("LogNormal"), rng) {} + LogNormalDistribution::LogNormalDistribution(const double tau, const double sigma) : BaseDistribution(std::string("LogNormal"), tau), m_sigma(sigma), m_distribution(std::log(tau), sigma) {} + LogNormalDistribution::LogNormalDistribution() : BaseDistribution(std::string("LogNormal")) {} void LogNormalDistribution::setParameters(const std::unordered_map ¶meters) { m_sigma = parameters.at("sigma"); @@ -16,13 +16,13 @@ namespace times { return parameter; } - void LogNormalDistribution::initialize() { + void LogNormalDistribution::initialize(std::mt19937_64& rng) { m_distribution = std::lognormal_distribution(std::log(m_tau), m_sigma); - m_tau_jump = m_distribution(m_rng); + m_tau_jump = m_distribution(rng); } - void LogNormalDistribution::draw_tau() { - m_tau_jump = m_distribution(m_rng); + std::unique_ptr LogNormalDistribution::clone() const { + return std::make_unique(*this); } std::string LogNormalDistribution::toString() const { diff --git a/src/times/lognormal.h b/src/times/lognormal.h index b2991ac..09f789b 100644 --- a/src/times/lognormal.h +++ b/src/times/lognormal.h @@ -7,16 +7,16 @@ namespace times { class LogNormalDistribution final : public BaseDistribution { public: - LogNormalDistribution(double, double, std::mt19937_64&); - explicit LogNormalDistribution(std::mt19937_64 &rng); + LogNormalDistribution(double, double); + LogNormalDistribution(); void setParameters(const std::unordered_map &) override; [[nodiscard]] std::unordered_map getParameters() const override; [[nodiscard]] std::string toString() const override; + [[nodiscard]] std::unique_ptr clone() const override; - void initialize() override; - void draw_tau() override; + void initialize(std::mt19937_64& rng) override; private: double m_sigma{1};