Use unique_ptr instead of raw pointer

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);
+    auto motion = motions::BaseMotion::createFromInput(args.motion_type);
-    times::BaseDistribution *dist = times::BaseDistribution::createFromInput(args.distribution_type, rng);
+    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;

src/motions/base.cpp

@@ -11,19 +11,16 @@
 #include "nsiteconejump.h"
 #include "sixsitejump.h"
 #include "rjoac.h"
+#include "conewobble.h"
 
 #include <stdexcept>
 
 
 
 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) {
+    BaseMotion::BaseMotion(std::string name, const double delta, const double eta) : m_name(std::move(name)), m_delta(delta), m_eta(eta) {}
-        m_uni_dist = std::uniform_real_distribution(0., 1.);
-    }
 
-    BaseMotion::BaseMotion(std::string name, std::mt19937_64& rng) : m_name(std::move(name)), m_rng(rng) {
+    BaseMotion::BaseMotion(std::string name) : m_name(std::move(name)) {}
-        m_uni_dist = std::uniform_real_distribution(0., 1.);
-    }
 
     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() {
+    coordinates::SphericalPos BaseMotion::draw_position(std::mt19937_64& rng) {
-        const double cos_theta = 1 - 2 * m_uni_dist(m_rng);
+        const double cos_theta = 1 - 2 * m_uni_dist(rng);
-        const double phi = 2.0 * M_PI * m_uni_dist(m_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> BaseMotion::createFromInput(const std::string& input) {
         if (input == "FourSiteTetrahedral")
-            return new FourSiteTetrahedron(rng);
+            return std::make_unique<FourSiteTetrahedron>();
 
         if (input == "SixSiteOctahedralC3")
-            return new SixSiteOctahedronC3(rng);
+            return std::make_unique<SixSiteOctahedronC3>();
 
         if (input == "IsotropicAngle")
-            return new SmallAngle(rng);
+            return std::make_unique<SmallAngle>();
 
         if (input == "RandomJump")
-            return new RandomJump(rng);
+            return std::make_unique<RandomJump>();
 
         if (input == "BimodalAngle")
-            return new BimodalAngle(rng);
+            return std::make_unique<BimodalAngle>();
 
         if (input == "NSiteConeJump")
-            return new NSiteJumpOnCone(rng);
+            return std::make_unique<NSiteJumpOnCone>();
 
         if (input == "RandomJumpOnCone")
-            return new RandomJumpOnCone(rng);
+            return std::make_unique<RandomJumpOnCone>();
+
+        if (input == "ConeWobble")
+            return std::make_unique<WobbleCone>();
 
         throw std::invalid_argument("Invalid input " + input);
     }

src/motions/base.h

@@ -3,6 +3,7 @@
 
 #include "coordinates.h"
 
+#include <memory>
 #include <random>
 #include <unordered_map>
 
@@ -11,15 +12,16 @@ namespace motions {
     public:
         virtual ~BaseMotion() = default;
 
-        BaseMotion(std::string, double, double, std::mt19937_64&);
+        BaseMotion(std::string, double, double);
-        explicit BaseMotion(std::string, std::mt19937_64&);
+        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 void initialize(std::mt19937_64& rng) = 0;
-        virtual double jump() = 0;
+        virtual double jump(std::mt19937_64& rng) = 0;
+        [[nodiscard]] virtual std::unique_ptr<BaseMotion> clone() const = 0;
 
         virtual void setParameters(const std::unordered_map<std::string, double>&);
         [[nodiscard]] virtual std::unordered_map<std::string, double> 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<BaseMotion> 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{0., 1.};
-        std::uniform_real_distribution<> m_uni_dist;
         double m_initial_omega{0.};
     };
 

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) :
+    BimodalAngle::BimodalAngle(const double delta, const double eta, const double angle1, const double angle2, const double prob) :
-        BaseMotion(std::string("BimodalAngle"), delta, eta, rng),
+        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() {
+    void BimodalAngle::initialize(std::mt19937_64& rng) {
-        m_prev_pos = draw_position();
+        m_prev_pos = draw_position(rng);
         m_initial_omega = omega_q(m_prev_pos);
     }
 
-    double BimodalAngle::jump() {
+    double BimodalAngle::jump(std::mt19937_64& rng) {
-        const double angle = m_uni_dist(m_rng) < m_prob ? m_angle1 : m_angle2;
+        const double angle = m_uni_dist(rng) < m_prob ? m_angle1 : m_angle2;
-        const double gamma{2 * M_PI * m_uni_dist(m_rng)};
+        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<BaseMotion> BimodalAngle::clone() const {
+        return std::make_unique<BimodalAngle>(*this);
+    }
+
     void BimodalAngle::setParameters(const std::unordered_map<std::string, double> &parameter) {
         BaseMotion::setParameters(parameter);
 

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& );
+        BimodalAngle(double, double, double, double, double);
-        explicit BimodalAngle(std::mt19937_64&);
+        BimodalAngle();
 
-        void initialize() override;
+        void initialize(std::mt19937_64& rng) override;
-        double jump() override;
+        double jump(std::mt19937_64& rng) override;
+        [[nodiscard]] std::unique_ptr<BaseMotion> clone() const override;
         void setParameters(const std::unordered_map<std::string, double> &) override;
         [[nodiscard]] std::unordered_map<std::string, double> getParameters() const override;
         [[nodiscard]] std::string toString() const override;

src/motions/conewobble.cpp

@@ -2,23 +2,24 @@
 #include "conewobble.h"
 #include "coordinates.h"
 
-#include <random>
-#include <string>
-
 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(const double delta, const double eta, const double chi) : BaseMotion("Wobble in Cone", delta, eta), m_angle(chi) {}
-    WobbleCone::WobbleCone(std::mt19937_64 &rng) : BaseMotion("Wobble in Cone", rng) {}
+    WobbleCone::WobbleCone() : BaseMotion("Wobble in Cone") {}
 
-    void WobbleCone::initialize() {
+    void WobbleCone::initialize(std::mt19937_64& rng) {
-        m_axis = draw_position();
+        m_axis = draw_position(rng);
     }
 
-    double WobbleCone::jump() {
+    double WobbleCone::jump(std::mt19937_64& rng) {
-        const double real_angle = m_uni_dist(m_rng) * m_angle;
+        const double real_angle = m_uni_dist(rng) * m_angle;
-        const double phi = 2 * M_PI * m_uni_dist(m_rng);
+        const double phi = 2 * M_PI * m_uni_dist(rng);
         return omega_q(rotate(m_axis, real_angle, phi));
     }
 
+    std::unique_ptr<BaseMotion> WobbleCone::clone() const {
+        return std::make_unique<WobbleCone>(*this);
+    }
+
     void WobbleCone::setParameters(const std::unordered_map<std::string, double> &parameters) {
         BaseMotion::setParameters(parameters);
         m_angle = parameters.at("angle");

src/motions/conewobble.h

@@ -4,22 +4,18 @@
 #include "base.h"
 #include "coordinates.h"
 
-#include <random>
-
 namespace motions {
     class WobbleCone final: public BaseMotion {
     public:
-        WobbleCone(double, double, double, std::mt19937_64&);
+        WobbleCone(double, double, double);
-        explicit WobbleCone(std::mt19937_64&);
+        WobbleCone();
 
-        void initialize() override;
+        void initialize(std::mt19937_64& rng) override;
-
+        double jump(std::mt19937_64& rng) override;
-        double jump() override;
+        [[nodiscard]] std::unique_ptr<BaseMotion> clone() const override;
 
         void setParameters(const std::unordered_map<std::string, double> &) override;
-
         [[nodiscard]] std::unordered_map<std::string, double> getParameters() const override;
-
         [[nodiscard]] std::string toString() const override;
 
     private:

src/motions/foursitejump.cpp

@@ -1,34 +1,36 @@
 #include "foursitejump.h"
 #include "coordinates.h"
 
-#include <random>
-
 namespace motions {
-    FourSiteTetrahedron::FourSiteTetrahedron(const double delta, const double eta, std::mt19937_64& rng) :
+    FourSiteTetrahedron::FourSiteTetrahedron(const double delta, const double eta) :
-        BaseMotion(std::string{"FourSiteTetrahedral"}, delta, eta, rng) {}
+        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() {
+    void FourSiteTetrahedron::initialize(std::mt19937_64& rng) {
-        const auto pos = draw_position();
+        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() {
+    double FourSiteTetrahedron::jump(std::mt19937_64& rng) {
-        m_corner_idx += m_chooser(m_rng);
+        m_corner_idx += m_chooser(rng);
         m_corner_idx %= 4;
 
         return m_corners[m_corner_idx];
     }
 
+    std::unique_ptr<BaseMotion> FourSiteTetrahedron::clone() const {
+        return std::make_unique<FourSiteTetrahedron>(*this);
+    }
+
     std::string FourSiteTetrahedron::toString() const {
         return {"FourSiteTetrahedral"};
     }

src/motions/foursitejump.h

@@ -3,18 +3,18 @@
 
 #include "base.h"
 
-#include <random>
 #include <cmath>
 #include <array>
 
 namespace motions {
     class FourSiteTetrahedron final : public BaseMotion {
     public:
-        FourSiteTetrahedron(double, double, std::mt19937_64&);
+        FourSiteTetrahedron(double, double);
-        explicit FourSiteTetrahedron(std::mt19937_64&);
+        FourSiteTetrahedron();
 
-        void initialize() override;
+        void initialize(std::mt19937_64& rng) override;
-        double jump() override;
+        double jump(std::mt19937_64& rng) override;
+        [[nodiscard]] std::unique_ptr<BaseMotion> clone() const override;
 
         [[nodiscard]] std::string toString() const override;
 

src/motions/isosmallangle.cpp

@@ -1,26 +1,28 @@
 #include "isosmallangle.h"
 #include "coordinates.h"
 
-#include <iostream>
-
 
 namespace motions {
-    SmallAngle::SmallAngle(const double delta, const double eta, const double chi, std::mt19937_64 &rng) :
+    SmallAngle::SmallAngle(const double delta, const double eta, const double chi) :
-        BaseMotion(std::string("IsotropicAngle"), delta, eta, rng), m_chi(chi * M_PI / 180.0) {}
+        BaseMotion(std::string("IsotropicAngle"), delta, eta), m_chi(chi * M_PI / 180.0) {}
-    SmallAngle::SmallAngle(std::mt19937_64 &rng) : BaseMotion(std::string("IsotropicAngle"), rng) {}
+    SmallAngle::SmallAngle() : BaseMotion(std::string("IsotropicAngle")) {}
 
-    void SmallAngle::initialize() {
+    void SmallAngle::initialize(std::mt19937_64& rng) {
-        m_prev_pos = draw_position();
+        m_prev_pos = draw_position(rng);
         m_initial_omega = omega_q(m_prev_pos);
     }
 
-    double SmallAngle::jump() {
+    double SmallAngle::jump(std::mt19937_64& rng) {
-        const double gamma{2 * M_PI * m_uni_dist(m_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<BaseMotion> SmallAngle::clone() const {
+        return std::make_unique<SmallAngle>(*this);
+    }
+
     void SmallAngle::setParameters(const std::unordered_map<std::string, double> &parameters) {
         m_chi = parameters.at("angle") * M_PI / 180.0;
         BaseMotion::setParameters(parameters);

src/motions/isosmallangle.h

@@ -8,11 +8,12 @@
 namespace motions {
     class SmallAngle final : public BaseMotion {
     public:
-        SmallAngle(double, double, double, std::mt19937_64& );
+        SmallAngle(double, double, double);
-        explicit SmallAngle(std::mt19937_64&);
+        SmallAngle();
 
-        void initialize() override;
+        void initialize(std::mt19937_64& rng) override;
-        double jump() override;
+        double jump(std::mt19937_64& rng) override;
+        [[nodiscard]] std::unique_ptr<BaseMotion> clone() const override;
         void setParameters(const std::unordered_map<std::string, double> &) override;
         [[nodiscard]] std::unordered_map<std::string, double> getParameters() const override;
         [[nodiscard]] std::string toString() const override;

src/motions/nsiteconejump.cpp

@@ -1,43 +1,42 @@
-//
-// Created by dominik on 12/14/24.
-//
 
 #include "nsiteconejump.h"
 #include "coordinates.h"
 
 #include <cmath>
-#include <iostream>
-#include <ostream>
 #include <vector>
 
 namespace motions {
-    NSiteJumpOnCone::NSiteJumpOnCone(const double delta, const double eta, const int num_sites, const double chi, std::mt19937_64 &rng) :
+    NSiteJumpOnCone::NSiteJumpOnCone(const double delta, const double eta, const int num_sites, const double chi) :
-    BaseMotion("NSiteJumpOnCone", delta, eta, rng),
+    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<double>(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() {
+    double NSiteJumpOnCone::jump(std::mt19937_64& rng) {
-        m_cone_idx += m_chooser(m_rng);
+        m_cone_idx += m_chooser(rng);
         m_cone_idx %= m_num_sites;
 
         return m_sites[m_cone_idx];
     }
 
+    std::unique_ptr<BaseMotion> NSiteJumpOnCone::clone() const {
+        return std::make_unique<NSiteJumpOnCone>(*this);
+    }
+
     void NSiteJumpOnCone::setParameters(const std::unordered_map<std::string, double> &parameters) {
         BaseMotion::setParameters(parameters);
         m_num_sites = static_cast<int>(parameters.at("num_sites"));

src/motions/nsiteconejump.h

@@ -3,17 +3,17 @@
 
 #include "base.h"
 
-#include <random>
 #include <vector>
 
 namespace motions {
     class NSiteJumpOnCone final : public BaseMotion {
     public:
-        NSiteJumpOnCone(double, double, int, double, std::mt19937_64&);
+        NSiteJumpOnCone(double, double, int, double);
-        explicit NSiteJumpOnCone(std::mt19937_64&);
+        NSiteJumpOnCone();
 
-        void initialize() override;
+        void initialize(std::mt19937_64& rng) override;
-        double jump() override;
+        double jump(std::mt19937_64& rng) override;
+        [[nodiscard]] std::unique_ptr<BaseMotion> clone() const override;
 
         [[nodiscard]] std::string toString() const override;
 

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() {
+    std::unique_ptr<BaseMotion> RandomJump::clone() const {
-        m_initial_omega = RandomJump::jump();
+        return std::make_unique<RandomJump>(*this);
     }
 
-    double RandomJump::jump() {
+    void RandomJump::initialize(std::mt19937_64& rng) {
-        return omega_q(draw_position());
+        m_initial_omega = RandomJump::jump(rng);
+    }
+
+    double RandomJump::jump(std::mt19937_64& rng) {
+        return omega_q(draw_position(rng));
     }
 }

src/motions/random.h

@@ -3,18 +3,18 @@
 
 
 #include "base.h"
-#include <random>
 
 namespace motions {
     class RandomJump final : public BaseMotion {
     public:
-        RandomJump(double, double, std::mt19937_64&);
+        RandomJump(double, double);
-        explicit RandomJump(std::mt19937_64&);
+        RandomJump();
 
         [[nodiscard]] std::string toString() const override;
+        [[nodiscard]] std::unique_ptr<BaseMotion> clone() const override;
 
-        void initialize() override;
+        void initialize(std::mt19937_64& rng) override;
-        double jump() override;
+        double jump(std::mt19937_64& rng) override;
     };
 }
 

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(const double delta, const double eta, const double chi) : BaseMotion("RJ on a Cone", delta, eta), m_angle(chi) {}
-    RandomJumpOnCone::RandomJumpOnCone(std::mt19937_64 &rng) : BaseMotion("RJ on a Cone", rng) {}
+    RandomJumpOnCone::RandomJumpOnCone() : BaseMotion("RJ on a Cone") {}
 
-    void RandomJumpOnCone::initialize() {
+    void RandomJumpOnCone::initialize(std::mt19937_64& rng) {
-        m_axis = draw_position();
+        m_axis = draw_position(rng);
     }
 
-    double RandomJumpOnCone::jump() {
+    double RandomJumpOnCone::jump(std::mt19937_64& rng) {
-        const double phi = 2 * M_PI * m_uni_dist(m_rng);
+        const double phi = 2 * M_PI * m_uni_dist(rng);
         return omega_q(rotate(m_axis, m_angle, phi));
     }
 
+    std::unique_ptr<BaseMotion> RandomJumpOnCone::clone() const {
+        return std::make_unique<RandomJumpOnCone>(*this);
+    }
+
     void RandomJumpOnCone::setParameters(const std::unordered_map<std::string, double> &parameters) {
         BaseMotion::setParameters(parameters);
         m_angle = parameters.at("angle");

src/motions/rjoac.h

@@ -4,22 +4,18 @@
 #include "base.h"
 #include "coordinates.h"
 
-#include <random>
-
 namespace motions {
     class RandomJumpOnCone final: public BaseMotion {
     public:
-        RandomJumpOnCone(double, double, double, std::mt19937_64&);
+        RandomJumpOnCone(double, double, double);
-        explicit RandomJumpOnCone(std::mt19937_64&);
+        RandomJumpOnCone();
 
-        void initialize() override;
+        void initialize(std::mt19937_64& rng) override;
-
+        double jump(std::mt19937_64& rng) override;
-        double jump() override;
+        [[nodiscard]] std::unique_ptr<BaseMotion> clone() const override;
 
         void setParameters(const std::unordered_map<std::string, double> &) override;
-
         [[nodiscard]] std::unordered_map<std::string, double> getParameters() const override;
-
         [[nodiscard]] std::string toString() const override;
 
     private:

src/motions/sixsitejump.cpp

@@ -1,21 +1,17 @@
 #include "sixsitejump.h"
 
-#include <iomanip>
-#include <iostream>
-
 #include "coordinates.h"
 
 namespace motions {
-    SixSiteOctahedronC3::SixSiteOctahedronC3(const double delta, const double eta, const double chi, std::mt19937_64& rng) :
+    SixSiteOctahedronC3::SixSiteOctahedronC3(const double delta, const double eta, const double chi) :
-        m_chi{chi*M_PI/180.},
+        BaseMotion(std::string{"SixSiteOctahedral"}, delta, eta),
-        BaseMotion(std::string{"SixSiteOctahedral"}, delta, eta, rng) {}
+        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() {
+    void SixSiteOctahedronC3::initialize(std::mt19937_64& rng) {
-        const coordinates::SphericalPos c3_axis = draw_position();
+        const coordinates::SphericalPos c3_axis = draw_position(rng);
-        const auto [x, y, z] = spherical_to_xyz(c3_axis);
+        const double alpha = 2. * M_PI * m_uni_dist(rng);
-        const double alpha = 2. * M_PI * m_uni_dist(m_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() {
+    double SixSiteOctahedronC3::jump(std::mt19937_64& rng) {
-        m_corner_idx += m_chooser(m_rng);
+        m_corner_idx += m_chooser(rng);
         m_corner_idx %= 6;
 
         return m_corners[m_corner_idx];
     }
 
+    std::unique_ptr<BaseMotion> SixSiteOctahedronC3::clone() const {
+        return std::make_unique<SixSiteOctahedronC3>(*this);
+    }
+
     std::string SixSiteOctahedronC3::toString() const {
         return {"SixSiteOctahedral/angle=" + std::to_string(m_chi / M_PI * 180.)};
     }

src/motions/sixsitejump.h

@@ -4,23 +4,23 @@
 
 #include "base.h"
 
-#include <random>
 #include <cmath>
 #include <array>
 
 namespace motions {
     class SixSiteOctahedronC3 final : public BaseMotion {
     public:
-        SixSiteOctahedronC3(double, double, double, std::mt19937_64&);
+        SixSiteOctahedronC3(double, double, double);
-        explicit SixSiteOctahedronC3(std::mt19937_64&);
+        SixSiteOctahedronC3();
 
-        void initialize() override;
+        void initialize(std::mt19937_64& rng) override;
-        double jump() override;
+        double jump(std::mt19937_64& rng) override;
+        [[nodiscard]] std::unique_ptr<BaseMotion> 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<double, 6> m_corners{};
         int m_corner_idx{0};

src/sims.cpp

@@ -19,7 +19,8 @@ void run_spectrum(
     std::unordered_map<std::string, double>& parameter,
     std::unordered_map<std::string, double> optional,
     motions::BaseMotion& motion,
-    times::BaseDistribution& dist
+    times::BaseDistribution& dist,
+    std::mt19937_64& rng
     ) {
     const int num_walker = static_cast<int>(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<double> traj_time{};
+        auto traj = make_trajectory(motion, dist, tmax, rng);
-        std::vector<double> traj_phase{};
-        std::vector<double> traj_omega{};
-
-        make_trajectory(motion, dist, tmax, traj_time, traj_phase, traj_omega);
 
         for (auto& [t_echo_j, fid_j] : fid_dict) {
             // get phase at echo pulse
-            int current_pos = nearest_index(traj_time, t_echo_j, 0);
+            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);
+            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);
+                current_pos = nearest_index(traj.time, real_time, current_pos);
-                const double phase_acq = lerp(traj_time, traj_phase, 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<std::string, double>& parameter,
     std::unordered_map<std::string, double> optional,
     motions::BaseMotion& motion,
-    times::BaseDistribution& dist
+    times::BaseDistribution& dist,
+    std::mt19937_64& rng
     ) {
     const int num_walker = static_cast<int>(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<double> traj_time{};
+        auto traj = make_trajectory(motion, dist, tmax, rng);
-        std::vector<double> traj_phase{};
-        std::vector<double> traj_omega{};
-
-        make_trajectory(motion, dist, tmax, traj_time, traj_phase, traj_omega);
 
         int f2_pos = 0;
         for (int f2_idx=0; f2_idx < num_mix_times; f2_idx++) {
             const double t_mix_f2 = mixing_times[f2_idx];
-            f2_pos = nearest_index(traj_time, t_mix_f2, f2_pos);
+            f2_pos = nearest_index(traj.time, t_mix_f2, f2_pos);
-            f2[f2_idx] += traj_omega[f2_pos] * motion.getInitOmega() / num_walker;
+            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);
+            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 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);
+                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);
+                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);
+                current_pos = nearest_index(traj.time, time_pulse4, current_pos);
-                const double phase_4pulse = lerp(traj_time, traj_phase, 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);
+                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;
+                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<double>& out_time,
+    std::mt19937_64& rng
-    std::vector<double>& out_phase,
-    std::vector<double>& out_omega
     ) {
     // Starting position
     double t_passed = 0;
     double phase = 0;
 
-    motion.initialize();
+    motion.initialize(rng);
-    dist.initialize();
+    dist.initialize(rng);
 
     double omega = motion.getInitOmega();
 
-    out_time.emplace_back(t_passed);
+    Trajectory traj;
-    out_phase.emplace_back(phase);
+    traj.time.emplace_back(t_passed);
-    out_omega.emplace_back(omega);
+    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);
+        traj.time.emplace_back(t_passed);
-        out_phase.emplace_back(phase);
+        traj.phase.emplace_back(phase);
-        out_omega.emplace_back(omega);
+        traj.omega.emplace_back(omega);
     }
+
+    return traj;
 }
 
 

src/sims.h

@@ -7,38 +7,19 @@
 #include <unordered_map>
 #include <string>
 #include <chrono>
+#include <vector>
 
-/**
+struct Trajectory {
- * @brief Run simulation for spectra
+    std::vector<double> time;
- *
+    std::vector<double> phase;
- * @param parameter Dictionary of parameter for simulation
+    std::vector<double> omega;
- * @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<std::string, double>& parameter, std::unordered_map<std::string, double> optional, motions::BaseMotion& motion, times::BaseDistribution& dist);
 
-/**
+void run_spectrum(std::unordered_map<std::string, double>& parameter, std::unordered_map<std::string, double> optional, motions::BaseMotion& motion, times::BaseDistribution& dist, std::mt19937_64& rng);
- * @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<std::string, double>& parameter, std::unordered_map<std::string, double> optional, motions::BaseMotion& motion, times::BaseDistribution& dist);
 
-/**
+void run_ste(std::unordered_map<std::string, double>& parameter, std::unordered_map<std::string, double> optional, motions::BaseMotion& motion, times::BaseDistribution& dist, std::mt19937_64& rng);
- * @brief Create trajectory of a single walker
+
- *
+Trajectory make_trajectory(motions::BaseMotion& motion, times::BaseDistribution& dist, double t_max, std::mt19937_64& rng);
- * @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<double>& out_time, std::vector<double>& out_phase, std::vector<double>& out_omega);
 
 std::chrono::system_clock::time_point printStart(std::unordered_map<std::string, double> &optional);
 void printEnd(std::chrono::system_clock::time_point start);

src/times/base.cpp

@@ -5,12 +5,12 @@
 #include <stdexcept>
 
 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 {
+    double BaseDistribution::tau_wait(std::mt19937_64& rng) const {
-        return std::exponential_distribution(1./m_tau_jump)(m_rng);
+        return std::exponential_distribution(1./m_tau_jump)(rng);
     }
 
     void BaseDistribution::setParameters(const std::unordered_map<std::string, double> &parameters) {
@@ -24,12 +24,12 @@ namespace times {
     }
 
 
-    BaseDistribution* BaseDistribution::createFromInput(const std::string& input, std::mt19937_64& rng) {
+    std::unique_ptr<BaseDistribution> BaseDistribution::createFromInput(const std::string& input) {
         if (input == "Delta")
-            return new DeltaDistribution(rng);
+            return std::make_unique<DeltaDistribution>();
 
         if (input == "LogNormal")
-            return new LogNormalDistribution(rng);
+            return std::make_unique<LogNormalDistribution>();
 
         throw std::invalid_argument("Invalid input " + input);
     }

src/times/base.h

@@ -1,6 +1,7 @@
 #ifndef RWSIM_TIMESBASE_H
 #define RWSIM_TIMESBASE_H
 
+#include <memory>
 #include <random>
 #include <unordered_map>
 
@@ -9,8 +10,8 @@ namespace times {
     public:
         virtual ~BaseDistribution() = default;
 
-        BaseDistribution(std::string, double, std::mt19937_64&);
+        BaseDistribution(std::string, double);
-        explicit BaseDistribution(std::string, std::mt19937_64&);
+        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<std::string, double>&);
         [[nodiscard]] virtual std::unordered_map<std::string, double> getParameters() const;
 
-        virtual void initialize() = 0;
+        virtual void initialize(std::mt19937_64& rng) = 0;
-        virtual void draw_tau() = 0;
+        [[nodiscard]] double tau_wait(std::mt19937_64& rng) const;
-        [[nodiscard]] double tau_wait() const;
+        [[nodiscard]] virtual std::unique_ptr<BaseDistribution> 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<BaseDistribution> 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;
     };
 }
 

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(const double tau) : BaseDistribution(std::string("Delta"), tau) {}
-    DeltaDistribution::DeltaDistribution(std::mt19937_64& rng) : BaseDistribution(std::string("Delta"), rng) {}
+    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<BaseDistribution> DeltaDistribution::clone() const {
+        return std::make_unique<DeltaDistribution>(*this);
+    }
 
     std::string DeltaDistribution::toString() const {
         return {"Delta/tau=" + std::to_string(m_tau)};

src/times/delta.h

@@ -6,11 +6,11 @@
 namespace times {
     class DeltaDistribution final : public BaseDistribution {
     public:
-        DeltaDistribution(double, std::mt19937_64&);
+        explicit DeltaDistribution(double);
-        explicit DeltaDistribution(std::mt19937_64 &rng);
+        DeltaDistribution();
 
-        void initialize() override;
+        void initialize(std::mt19937_64& rng) override;
-        void draw_tau() override;
+        [[nodiscard]] std::unique_ptr<BaseDistribution> clone() const override;
 
         [[nodiscard]] std::string toString() const override;
     };

src/times/lognormal.cpp

@@ -2,8 +2,8 @@
 #include <cmath>
 
 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(const double tau, const double sigma) : BaseDistribution(std::string("LogNormal"), tau), m_sigma(sigma), m_distribution(std::log(tau), sigma) {}
-    LogNormalDistribution::LogNormalDistribution(std::mt19937_64& rng) : BaseDistribution(std::string("LogNormal"), rng) {}
+    LogNormalDistribution::LogNormalDistribution() : BaseDistribution(std::string("LogNormal")) {}
 
     void LogNormalDistribution::setParameters(const std::unordered_map<std::string, double> &parameters) {
         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() {
+    std::unique_ptr<BaseDistribution> LogNormalDistribution::clone() const {
-        m_tau_jump = m_distribution(m_rng);
+        return std::make_unique<LogNormalDistribution>(*this);
     }
 
     std::string LogNormalDistribution::toString() const {

src/times/lognormal.h

@@ -7,16 +7,16 @@
 namespace times {
     class LogNormalDistribution final : public BaseDistribution {
     public:
-        LogNormalDistribution(double, double, std::mt19937_64&);
+        LogNormalDistribution(double, double);
-        explicit LogNormalDistribution(std::mt19937_64 &rng);
+        LogNormalDistribution();
 
         void setParameters(const std::unordered_map<std::string, double> &) override;
         [[nodiscard]] std::unordered_map<std::string, double> getParameters() const override;
 
         [[nodiscard]] std::string toString() const override;
+        [[nodiscard]] std::unique_ptr<BaseDistribution> clone() const override;
 
-        void initialize() override;
+        void initialize(std::mt19937_64& rng) override;
-        void draw_tau() override;
 
     private:
         double m_sigma{1};