add namespaces and cmakelists for sub-dirs

This commit is contained in:
Dominik Demuth 2024-12-13 14:03:33 +01:00
parent e90c4c9543
commit e6331749b2
32 changed files with 604 additions and 440 deletions

View File

@ -1,35 +1,17 @@
cmake_minimum_required(VERSION 3.18)
set(CMAKE_CXX_STANDARD 17) add_subdirectory(times)
add_subdirectory(motions)
add_subdirectory(utils)
add_library(simulation STATIC sims.cpp sims.h)
target_link_libraries(simulation PRIVATE utils)
add_executable(rwsim main.cpp add_executable(
utils/functions.h rwsim
utils/functions.cpp main.cpp
utils/io.cpp
utils/io.h
motions/base.cpp
motions/base.h
motions/random.cpp
motions/random.h
times/base.cpp
times/base.h
times/delta.cpp
times/delta.h
simulation/sims.cpp
simulation/sims.h
utils/ranges.cpp
utils/ranges.h
motions/tetrahedral.cpp
motions/tetrahedral.h
motions/isosmallangle.cpp
motions/isosmallangle.h
motions/coordinates.cpp
motions/coordinates.h
motions/bimodalangle.cpp
motions/bimodalangle.h
times/lognormal.cpp
times/lognormal.h
) )
target_link_libraries(rwsim PUBLIC RWMotion RWTime utils simulation)
target_compile_options(rwsim PUBLIC -Werror -Wall -Wextra -Wconversion -O2) target_compile_options(rwsim PUBLIC -Werror -Wall -Wextra -Wconversion -O2)

View File

@ -1,10 +1,9 @@
#include "sims.h"
#include "utils/io.h" #include "utils/io.h"
#include "simulation/sims.h"
#include "motions/base.h" #include "motions/base.h"
#include "times/base.h" #include "times/base.h"
#include <iostream> #include <iostream>
#include <unordered_map> #include <unordered_map>
#include <random> #include <random>
@ -37,8 +36,8 @@ int main (const int argc, char *argv[]) {
std::random_device rd; std::random_device rd;
std::mt19937_64 rng(rd()); std::mt19937_64 rng(rd());
Motion *motion = Motion::createFromInput(args.motion_type, rng); motions::BaseMotion *motion = motions::BaseMotion::createFromInput(args.motion_type, rng);
Distribution *dist = Distribution::createFromInput(args.distribution_type, rng); times::BaseDistribution *dist = times::BaseDistribution::createFromInput(args.distribution_type, rng);
if (args.spectrum) { if (args.spectrum) {
run_spectrum(parameter, args.optional, *motion, *dist); run_spectrum(parameter, args.optional, *motion, *dist);
} }

View File

@ -0,0 +1,23 @@
# Create a library target for motions
add_library(
RWMotion STATIC
conewobble.cpp
conewobble.h
coordinates.cpp
coordinates.h
base.cpp
base.h
random.cpp
random.h
isosmallangle.cpp
isosmallangle.h
foursitejump.cpp
foursitejump.h
rjoac.cpp
rjoac.h
bimodalangle.cpp
bimodalangle.h
sixsitejump.cpp
sixsitejump.h
)

View File

@ -4,42 +4,47 @@
#include "bimodalangle.h" #include "bimodalangle.h"
#include "isosmallangle.h" #include "isosmallangle.h"
#include "random.h" #include "random.h"
#include "tetrahedral.h" #include "foursitejump.h"
#include <stdexcept> #include <stdexcept>
#include <map>
Motion::Motion(std::string name, const double delta, const double eta, std::mt19937_64& rng) : m_name(std::move(name)), m_delta(delta), m_eta(eta), m_rng(rng) { #include "sixsitejump.h"
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.); m_uni_dist = std::uniform_real_distribution(0., 1.);
} }
Motion::Motion(std::string name, std::mt19937_64& rng) : m_name(std::move(name)), m_rng(rng) { 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.); m_uni_dist = std::uniform_real_distribution(0., 1.);
} }
double Motion::omega_q(const double cos_theta, const double phi) const { double BaseMotion::omega_q(const double cos_theta, const double phi) const {
const double cos_theta_square = cos_theta * cos_theta; const double cos_theta_square = cos_theta * cos_theta;
const double sin_theta_square = 1. - cos_theta_square; const double sin_theta_square = 1. - cos_theta_square;
return M_PI * m_delta * (3. * cos_theta_square - 1. - m_eta * sin_theta_square * std::cos(2.*phi)); return M_PI * m_delta * (3. * cos_theta_square - 1. - m_eta * sin_theta_square * std::cos(2.*phi));
} }
double Motion::omega_q(const SphericalPos& pos) const { double BaseMotion::omega_q(const coordinates::SphericalPos& pos) const {
auto [cos_theta, phi] = pos; auto [cos_theta, phi] = pos;
return omega_q(cos_theta, phi); return omega_q(cos_theta, phi);
} }
SphericalPos Motion::draw_position() { coordinates::SphericalPos BaseMotion::draw_position() {
const double cos_theta = 1 - 2 * m_uni_dist(m_rng); const double cos_theta = 1 - 2 * m_uni_dist(m_rng);
const double phi = 2.0 * M_PI * m_uni_dist(m_rng); const double phi = 2.0 * M_PI * m_uni_dist(m_rng);
return {cos_theta, phi}; return {cos_theta, phi};
} }
Motion* Motion::createFromInput(const std::string& input, std::mt19937_64& rng) { BaseMotion* BaseMotion::createFromInput(const std::string& input, std::mt19937_64& rng) {
if (input == "TetrahedralJump") if (input == "FourSiteTetrahedral")
return new TetrahedralJump(rng); return new FourSiteTetrahedron(rng);
if (input == "SixSiteOctahedral")
return new SixSiteOctahedron(rng);
if (input == "IsotropicAngle") if (input == "IsotropicAngle")
return new SmallAngle(rng); return new SmallAngle(rng);
@ -53,19 +58,20 @@ Motion* Motion::createFromInput(const std::string& input, std::mt19937_64& rng)
throw std::invalid_argument("Invalid input " + input); throw std::invalid_argument("Invalid input " + input);
} }
void Motion::setParameters(const std::unordered_map<std::string, double> &parameters) { void BaseMotion::setParameters(const std::unordered_map<std::string, double> &parameters) {
m_delta = parameters.at("delta"); m_delta = parameters.at("delta");
m_eta = parameters.at("eta"); m_eta = parameters.at("eta");
} }
std::unordered_map<std::string, double> Motion::getParameters() const { std::unordered_map<std::string, double> BaseMotion::getParameters() const {
return std::unordered_map<std::string, double>{ return std::unordered_map<std::string, double>{
{"delta", m_delta}, {"delta", m_delta},
{"eta", m_eta} {"eta", m_eta}
}; };
} }
std::ostream& operator<<(std::ostream& os, const Motion& m) { std::ostream& operator<<(std::ostream& os, const BaseMotion& m) {
os << m.getName(); os << m.getName();
return os; return os;
} }
}

View File

@ -6,16 +6,17 @@
#include <random> #include <random>
#include <unordered_map> #include <unordered_map>
class Motion { namespace motions {
class BaseMotion {
public: public:
virtual ~Motion() = default; virtual ~BaseMotion() = default;
Motion(std::string, double, double, std::mt19937_64&); BaseMotion(std::string, double, double, std::mt19937_64&);
explicit Motion(std::string, std::mt19937_64&); explicit BaseMotion(std::string, std::mt19937_64&);
SphericalPos draw_position(); coordinates::SphericalPos draw_position();
[[nodiscard]] double omega_q(double, double) const; [[nodiscard]] double omega_q(double, double) const;
[[nodiscard]] double omega_q(const SphericalPos&) const; [[nodiscard]] double omega_q(const coordinates::SphericalPos&) const;
virtual void initialize() = 0; virtual void initialize() = 0;
virtual double jump() = 0; virtual double jump() = 0;
@ -31,7 +32,7 @@ public:
[[nodiscard]] virtual std::string toString() const = 0; [[nodiscard]] virtual std::string toString() const = 0;
static Motion* createFromInput(const std::string& input, std::mt19937_64& rng); static BaseMotion* createFromInput(const std::string& input, std::mt19937_64& rng);
protected: protected:
std::string m_name{"BaseMotion"}; std::string m_name{"BaseMotion"};
@ -42,6 +43,6 @@ protected:
double m_initial_omega{0.}; double m_initial_omega{0.};
}; };
std::ostream& operator<<(std::ostream& os, const Motion& m); std::ostream& operator<<(std::ostream& os, const BaseMotion& m);
}
#endif //RWSIM_MOTIONBASE_H #endif //RWSIM_MOTIONBASE_H

View File

@ -1,13 +1,15 @@
#include "bimodalangle.h" #include "bimodalangle.h"
#include "base.h" #include "base.h"
#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, std::mt19937_64 &rng) :
Motion(std::string("BimodalAngle"), delta, eta, rng), BaseMotion(std::string("BimodalAngle"), delta, eta, rng),
m_angle1(angle1 * M_PI / 180.0), m_angle1(angle1 * M_PI / 180.0),
m_angle2(angle2 * M_PI / 180.0), m_angle2(angle2 * M_PI / 180.0),
m_prob(prob) {}; m_prob(prob) {};
BimodalAngle::BimodalAngle(std::mt19937_64 &rng) : Motion(std::string("BimodalAngle"), rng) {} BimodalAngle::BimodalAngle(std::mt19937_64 &rng) : BaseMotion(std::string("BimodalAngle"), rng) {}
void BimodalAngle::initialize() { void BimodalAngle::initialize() {
m_prev_pos = draw_position(); m_prev_pos = draw_position();
@ -17,13 +19,13 @@ void BimodalAngle::initialize() {
double BimodalAngle::jump() { double BimodalAngle::jump() {
const double angle = m_uni_dist(m_rng) < m_prob ? m_angle1 : m_angle2; 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)}; const double gamma{2 * M_PI * m_uni_dist(m_rng)};
m_prev_pos = rotate(m_prev_pos, angle, gamma); m_prev_pos = coordinates::rotate(m_prev_pos, angle, gamma);
return omega_q(m_prev_pos); return omega_q(m_prev_pos);
} }
void BimodalAngle::setParameters(const std::unordered_map<std::string, double> &parameter) { void BimodalAngle::setParameters(const std::unordered_map<std::string, double> &parameter) {
Motion::setParameters(parameter); BaseMotion::setParameters(parameter);
m_angle1 = parameter.at("angle1") * M_PI / 180.; m_angle1 = parameter.at("angle1") * M_PI / 180.;
m_angle2 = parameter.at("angle2") * M_PI / 180.; m_angle2 = parameter.at("angle2") * M_PI / 180.;
@ -31,7 +33,7 @@ void BimodalAngle::setParameters(const std::unordered_map<std::string, double> &
} }
std::unordered_map<std::string, double> BimodalAngle::getParameters() const { std::unordered_map<std::string, double> BimodalAngle::getParameters() const {
auto parameter = Motion::getParameters(); auto parameter = BaseMotion::getParameters();
parameter["angle1"] = m_angle1 * 180 / M_PI; parameter["angle1"] = m_angle1 * 180 / M_PI;
parameter["angle2"] = m_angle2 * 180 / M_PI; parameter["angle2"] = m_angle2 * 180 / M_PI;
parameter["probality1"] = m_prob; parameter["probality1"] = m_prob;
@ -41,3 +43,4 @@ std::unordered_map<std::string, double> BimodalAngle::getParameters() const {
std::string BimodalAngle::toString() const { std::string BimodalAngle::toString() const {
return std::string{"BimodalAngle/angle1=" + std::to_string(m_angle1 * 180 / M_PI) + "/angle2=" + std::to_string(m_angle2 * 180 / M_PI) + "/probability1=" + std::to_string(m_prob)}; return std::string{"BimodalAngle/angle1=" + std::to_string(m_angle1 * 180 / M_PI) + "/angle2=" + std::to_string(m_angle2 * 180 / M_PI) + "/probability1=" + std::to_string(m_prob)};
} }
}

View File

@ -3,8 +3,10 @@
#define BIMODALANGLE_H #define BIMODALANGLE_H
#include "base.h" #include "base.h"
#include "coordinates.h"
class BimodalAngle : public Motion { namespace motions {
class BimodalAngle : public BaseMotion {
public: public:
BimodalAngle(double, double, double, double, double, std::mt19937_64& ); BimodalAngle(double, double, double, double, double, std::mt19937_64& );
explicit BimodalAngle(std::mt19937_64&); explicit BimodalAngle(std::mt19937_64&);
@ -19,7 +21,7 @@ protected:
double m_angle1{0}; double m_angle1{0};
double m_angle2{0}; double m_angle2{0};
double m_prob{0}; double m_prob{0};
SphericalPos m_prev_pos{0., 0.}; coordinates::SphericalPos m_prev_pos{0., 0.};
}; };
}
#endif //BIMODALANGLE_H #endif //BIMODALANGLE_H

View File

@ -3,6 +3,7 @@
#include <cmath> #include <cmath>
#include <iostream> #include <iostream>
namespace coordinates {
SphericalPos rotate(const SphericalPos& old_pos, const double alpha, const double beta) { SphericalPos rotate(const SphericalPos& old_pos, const double alpha, const double beta) {
const double sin_alpha{std::sin(alpha)}; const double sin_alpha{std::sin(alpha)};
const double cos_alpha{std::cos(alpha)}; const double cos_alpha{std::cos(alpha)};
@ -37,3 +38,4 @@ CartesianPos spherical_to_xyz(const SphericalPos& pos) {
SphericalPos xyz_to_spherical(const CartesianPos& pos) { SphericalPos xyz_to_spherical(const CartesianPos& pos) {
return {pos.z, std::atan2(pos.y, pos.x)}; return {pos.z, std::atan2(pos.y, pos.x)};
} }
}

View File

@ -2,6 +2,7 @@
#ifndef COORDINATES_H #ifndef COORDINATES_H
#define COORDINATES_H #define COORDINATES_H
namespace coordinates {
struct SphericalPos { struct SphericalPos {
double cos_theta; double cos_theta;
double phi; double phi;
@ -16,5 +17,6 @@ struct CartesianPos {
SphericalPos rotate(const SphericalPos&, double, double); SphericalPos rotate(const SphericalPos&, double, double);
CartesianPos spherical_to_xyz(const SphericalPos&); CartesianPos spherical_to_xyz(const SphericalPos&);
SphericalPos xyz_to_spherical(const CartesianPos&); SphericalPos xyz_to_spherical(const CartesianPos&);
}
#endif //COORDINATES_H #endif //COORDINATES_H

View File

@ -0,0 +1,35 @@
#include "foursitejump.h"
#include "coordinates.h"
#include <random>
namespace motions {
FourSiteTetrahedron::FourSiteTetrahedron(const double delta, const double eta, std::mt19937_64& rng) :
BaseMotion(std::string{"FourSiteTetrahedral"}, delta, eta, rng) {}
FourSiteTetrahedron::FourSiteTetrahedron(std::mt19937_64& rng) : BaseMotion(std::string{"FourSiteTetrahedral"}, rng) {}
void FourSiteTetrahedron::initialize() {
const auto pos = draw_position();
m_corners[0] = omega_q(pos);
const double alpha = 2. * M_PI * m_uni_dist(m_rng);
for (int i = 1; i<4; i++) {
auto corner_pos = coordinates::rotate(pos, m_beta, alpha + (i-1) * 2*M_PI/3.);
m_corners[i] = omega_q(corner_pos);
}
m_initial_omega = FourSiteTetrahedron::jump();
}
double FourSiteTetrahedron::jump() {
m_corner_idx += m_chooser(m_rng);
m_corner_idx %= 4;
return m_corners[m_corner_idx];
}
std::string FourSiteTetrahedron::toString() const {
return {"FourSiteTetrahedral"};
}
}

View File

@ -0,0 +1,31 @@
#ifndef RWSIM_MOTIONTETRAHEDRAL_H
#define RWSIM_MOTIONTETRAHEDRAL_H
#include "base.h"
#include <random>
#include <cmath>
#include <array>
namespace motions {
class FourSiteTetrahedron final : public BaseMotion {
public:
FourSiteTetrahedron(double, double, std::mt19937_64&);
explicit FourSiteTetrahedron(std::mt19937_64&);
void initialize() override;
double jump() override;
[[nodiscard]] std::string toString() const override;
private:
const double m_beta{std::acos(-1/3.)};
std::array<double, 4> m_corners{};
int m_corner_idx{0};
std::uniform_int_distribution<> m_chooser{1, 3};
};
}
#endif //RWSIM_MOTIONTETRAHEDRAL_H

View File

@ -4,10 +4,10 @@
#include <iostream> #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, std::mt19937_64 &rng) :
Motion(std::string("IsotropicAngle"), delta, eta, rng), m_chi(chi * M_PI / 180.0) {}; BaseMotion(std::string("IsotropicAngle"), delta, eta, rng), m_chi(chi * M_PI / 180.0) {};
SmallAngle::SmallAngle(std::mt19937_64 &rng) : Motion(std::string("IsotropicAngle"), rng) {} SmallAngle::SmallAngle(std::mt19937_64 &rng) : BaseMotion(std::string("IsotropicAngle"), rng) {}
void SmallAngle::initialize() { void SmallAngle::initialize() {
m_prev_pos = draw_position(); m_prev_pos = draw_position();
@ -16,18 +16,18 @@ void SmallAngle::initialize() {
double SmallAngle::jump() { double SmallAngle::jump() {
const double gamma{2 * M_PI * m_uni_dist(m_rng)}; const double gamma{2 * M_PI * m_uni_dist(m_rng)};
m_prev_pos = rotate(m_prev_pos, m_chi, gamma); m_prev_pos = coordinates::rotate(m_prev_pos, m_chi, gamma);
return omega_q(m_prev_pos); return omega_q(m_prev_pos);
} }
void SmallAngle::setParameters(const std::unordered_map<std::string, double> &parameters) { void SmallAngle::setParameters(const std::unordered_map<std::string, double> &parameters) {
m_chi = parameters.at("angle") * M_PI / 180.0; m_chi = parameters.at("angle") * M_PI / 180.0;
Motion::setParameters(parameters); BaseMotion::setParameters(parameters);
} }
std::unordered_map<std::string, double> SmallAngle::getParameters() const { std::unordered_map<std::string, double> SmallAngle::getParameters() const {
auto parameter = Motion::getParameters(); auto parameter = BaseMotion::getParameters();
parameter["angle"] = m_chi * 180 / M_PI; parameter["angle"] = m_chi * 180 / M_PI;
return parameter; return parameter;
} }
@ -35,3 +35,4 @@ std::unordered_map<std::string, double> SmallAngle::getParameters() const {
std::string SmallAngle::toString() const { std::string SmallAngle::toString() const {
return std::string{"IsotropicAngle/angle=" + std::to_string(m_chi * 180 / M_PI)}; return std::string{"IsotropicAngle/angle=" + std::to_string(m_chi * 180 / M_PI)};
} }
}

View File

@ -5,7 +5,8 @@
#include "base.h" #include "base.h"
#include "coordinates.h" #include "coordinates.h"
class SmallAngle final : public Motion { namespace motions {
class SmallAngle final : public BaseMotion {
public: public:
SmallAngle(double, double, double, std::mt19937_64& ); SmallAngle(double, double, double, std::mt19937_64& );
explicit SmallAngle(std::mt19937_64&); explicit SmallAngle(std::mt19937_64&);
@ -18,7 +19,7 @@ public:
private: private:
double m_chi{0}; double m_chi{0};
SphericalPos m_prev_pos{0., 0.}; coordinates::SphericalPos m_prev_pos{0., 0.};
}; };
}
#endif //RWSIM_MOTIONISOSMALLANGLE_H #endif //RWSIM_MOTIONISOSMALLANGLE_H

View File

@ -1,10 +1,10 @@
#include "random.h" #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, std::mt19937_64 &rng) : Motion(std::string("RandomJump"), delta, eta, rng) {} RandomJump::RandomJump(std::mt19937_64 &rng) : BaseMotion(std::string("RandomJump"), rng) {}
RandomJump::RandomJump(std::mt19937_64 &rng) : Motion(std::string("RandomJump"), rng) {}
std::string RandomJump::toString() const { std::string RandomJump::toString() const {
return {"RandomJump"}; return {"RandomJump"};
@ -17,3 +17,4 @@ void RandomJump::initialize() {
double RandomJump::jump() { double RandomJump::jump() {
return omega_q(draw_position()); return omega_q(draw_position());
} }
}

View File

@ -5,7 +5,8 @@
#include "base.h" #include "base.h"
#include <random> #include <random>
class RandomJump final : public Motion { namespace motions {
class RandomJump final : public BaseMotion {
public: public:
RandomJump(double, double, std::mt19937_64&); RandomJump(double, double, std::mt19937_64&);
explicit RandomJump(std::mt19937_64&); explicit RandomJump(std::mt19937_64&);
@ -15,5 +16,6 @@ public:
void initialize() override; void initialize() override;
double jump() override; double jump() override;
}; };
}
#endif //RWSIM_MOTIONRANDOMJUMP_H #endif //RWSIM_MOTIONRANDOMJUMP_H

32
src/motions/rjoac.cpp Normal file
View File

@ -0,0 +1,32 @@
#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) {}
void RandomJumpOnCone::initialize() {
m_axis = draw_position();
}
double RandomJumpOnCone::jump() {
const double phi = 2 * M_PI * m_uni_dist(m_rng);
return omega_q(coordinates::rotate(m_axis, m_angle, phi));
}
void RandomJumpOnCone::setParameters(const std::unordered_map<std::string, double> &parameters) {
BaseMotion::setParameters(parameters);
m_angle = parameters.at("angle");
}
std::unordered_map<std::string, double> RandomJumpOnCone::getParameters() const {
auto parameter = BaseMotion::getParameters();
parameter["angle"] = m_angle;
return parameter;
}
std::string RandomJumpOnCone::toString() const {
return std::string("RandomJumpOnCone/angle=") + std::to_string(m_angle);
}
}

31
src/motions/rjoac.h Normal file
View File

@ -0,0 +1,31 @@
#ifndef RJOAC_H
#define RJOAC_H
#include "base.h"
#include "coordinates.h"
#include <random>
namespace motions {
class RandomJumpOnCone: public BaseMotion {
public:
RandomJumpOnCone(double, double, double, std::mt19937_64&);
explicit RandomJumpOnCone(std::mt19937_64&);
void initialize() override;
double jump() 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:
double m_angle{0};
coordinates::SphericalPos m_axis{1, 0};
};
}
#endif //RJOAC_H

View File

@ -0,0 +1,34 @@
#include "sixsitejump.h"
#include "coordinates.h"
namespace motions {
SixSiteOctahedron::SixSiteOctahedron(const double delta, const double eta, std::mt19937_64& rng) :
BaseMotion(std::string{"SixSiteOctahedral"}, delta, eta, rng) {}
SixSiteOctahedron::SixSiteOctahedron(std::mt19937_64& rng) : BaseMotion(std::string{"SixSiteOctahedral"}, rng) {}
void SixSiteOctahedron::initialize() {
const auto pos = draw_position();
m_corners[0] = omega_q(pos);
m_corners[1] = omega_q(rotate(pos, M_PI, 0));
const double alpha = 2. * M_PI * m_uni_dist(m_rng);
for (int i = 2; i<6; i++) {
auto corner_pos = coordinates::rotate(pos, M_PI_2, alpha + (i-2) * M_PI_2);
m_corners[i] = omega_q(corner_pos);
}
m_initial_omega = SixSiteOctahedron::jump();
}
double SixSiteOctahedron::jump() {
m_corner_idx += m_chooser(m_rng);
m_corner_idx %= 6;
return m_corners[m_corner_idx];
}
std::string SixSiteOctahedron::toString() const {
return {"SixSiteOctahedral"};
}
}

32
src/motions/sixsitejump.h Normal file
View File

@ -0,0 +1,32 @@
#ifndef SIXSITEJUMP_H
#define SIXSITEJUMP_H
#include "base.h"
#include <random>
#include <cmath>
#include <array>
namespace motions {
class SixSiteOctahedron final : public BaseMotion {
public:
SixSiteOctahedron(double, double, std::mt19937_64&);
explicit SixSiteOctahedron(std::mt19937_64&);
void initialize() override;
double jump() override;
[[nodiscard]] std::string toString() const override;
private:
const double m_beta{std::acos(-1/3.)};
std::array<double, 6> m_corners{};
int m_corner_idx{0};
std::uniform_int_distribution<> m_chooser{1, 5};
};
}
#endif //SIXSITEJUMP_H

View File

@ -1,34 +0,0 @@
#include "tetrahedral.h"
#include <random>
#include "tetrahedral.h"
TetrahedralJump::TetrahedralJump(const double delta, const double eta, std::mt19937_64& rng) :
Motion(std::string{"FourSiteTetrahedral"}, delta, eta, rng) {}
TetrahedralJump::TetrahedralJump(std::mt19937_64& rng) : Motion(std::string{"FourSiteTetrahedral"}, rng) {}
void TetrahedralJump::initialize() {
const auto pos = draw_position();
m_corners[0] = omega_q(pos);
const double alpha = 2. * M_PI * m_uni_dist(m_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 = TetrahedralJump::jump();
}
double TetrahedralJump::jump() {
m_corner_idx += m_chooser(m_rng);
m_corner_idx %= 4;
return m_corners[m_corner_idx];
}
std::string TetrahedralJump::toString() const {
return {"FourSiteTetrahedral"};
}

View File

@ -1,29 +0,0 @@
#ifndef RWSIM_MOTIONTETRAHEDRAL_H
#define RWSIM_MOTIONTETRAHEDRAL_H
#include "base.h"
#include <random>
#include <cmath>
#include <array>
class TetrahedralJump final : public Motion {
public:
TetrahedralJump(double, double, std::mt19937_64&);
explicit TetrahedralJump(std::mt19937_64&);
void initialize() override;
double jump() override;
[[nodiscard]] std::string toString() const override;
private:
const double m_beta{std::acos(-1/3.)};
std::array<double, 4> m_corners{};
int m_corner_idx{0};
std::uniform_int_distribution<> m_chooser{1, 3};
};
#endif //RWSIM_MOTIONTETRAHEDRAL_H

View File

@ -1,9 +1,8 @@
#include "sims.h" #include "simulation/sims.h"
#include "../motions/base.h" #include "times/base.h"
#include "../times/base.h" #include "utils/functions.h"
#include "../utils/functions.h" #include "utils/ranges.h"
#include "../utils/ranges.h" #include "utils/io.h"
#include "../utils/io.h"
#include <iostream> #include <iostream>
#include <algorithm> #include <algorithm>
@ -15,11 +14,12 @@
#include <chrono> #include <chrono>
void run_spectrum( void run_spectrum(
std::unordered_map<std::string, double>& parameter, std::unordered_map<std::string, double>& parameter,
std::unordered_map<std::string, double> optional, std::unordered_map<std::string, double> optional,
Motion& motion, motions::BaseMotion& motion,
Distribution& dist times::BaseDistribution& dist
) { ) {
const int num_walker = static_cast<int>(parameter["num_walker"]); const int num_walker = static_cast<int>(parameter["num_walker"]);
@ -83,8 +83,8 @@ void run_spectrum(
void run_ste( void run_ste(
std::unordered_map<std::string, double>& parameter, std::unordered_map<std::string, double>& parameter,
std::unordered_map<std::string, double> optional, std::unordered_map<std::string, double> optional,
Motion& motion, motions::BaseMotion& motion,
Distribution& dist times::BaseDistribution& dist
) { ) {
const int num_walker = static_cast<int>(parameter[std::string("num_walker")]); const int num_walker = static_cast<int>(parameter[std::string("num_walker")]);
@ -175,8 +175,8 @@ void run_ste(
void make_trajectory( void make_trajectory(
Motion& motion, motions::BaseMotion& motion,
Distribution& dist, times::BaseDistribution& dist,
const double t_max, const double t_max,
std::vector<double>& out_time, std::vector<double>& out_time,
std::vector<double>& out_phase, std::vector<double>& out_phase,

View File

@ -1,8 +1,8 @@
#ifndef RWSIM_SIMS_H #ifndef RWSIM_SIMS_H
#define RWSIM_SIMS_H #define RWSIM_SIMS_H
#include "../motions/base.h" #include "motions/base.h"
#include "../times/base.h" #include "times/base.h"
#include <unordered_map> #include <unordered_map>
#include <string> #include <string>
@ -16,7 +16,7 @@
* @param motion Motion model * @param motion Motion model
* @param dist Distribution of correlation times * @param dist Distribution of correlation times
*/ */
void run_spectrum(std::unordered_map<std::string, double>& parameter, std::unordered_map<std::string, double> optional, Motion& motion, Distribution& dist); void run_spectrum(std::unordered_map<std::string, double>& parameter, std::unordered_map<std::string, double> optional, motions::BaseMotion& motion, times::BaseDistribution& dist);
/** /**
* @brief Run simulation for stimulated echoes * @brief Run simulation for stimulated echoes
@ -26,7 +26,7 @@ void run_spectrum(std::unordered_map<std::string, double>& parameter, std::unord
* @param motion Motion model * @param motion Motion model
* @param dist Distribution of correlation times * @param dist Distribution of correlation times
*/ */
void run_ste(std::unordered_map<std::string, double>& parameter, std::unordered_map<std::string, double> optional, Motion& motion, Distribution& dist); void run_ste(std::unordered_map<std::string, double>& parameter, std::unordered_map<std::string, double> optional, motions::BaseMotion& motion, times::BaseDistribution& dist);
/** /**
* @brief Create trajectory of a single walker * @brief Create trajectory of a single walker
@ -38,7 +38,7 @@ void run_ste(std::unordered_map<std::string, double>& parameter, std::unordered_
* @param out_phase Vector of phase between waiting times * @param out_phase Vector of phase between waiting times
* @param out_omega Vector of omega at jump time * @param out_omega Vector of omega at jump time
*/ */
void make_trajectory(Motion& motion, Distribution& dist, double t_max, std::vector<double>& out_time, std::vector<double>& out_phase, std::vector<double>& out_omega); 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); std::chrono::system_clock::time_point printStart(std::unordered_map<std::string, double> &optional);
void printEnd(std::chrono::system_clock::time_point start); void printEnd(std::chrono::system_clock::time_point start);

View File

@ -4,26 +4,27 @@
#include <stdexcept> #include <stdexcept>
Distribution::Distribution(std::string name, const double tau, std::mt19937_64 &rng) : m_name(std::move(name)), m_tau(tau), m_tau_jump(tau), m_rng(rng) {} 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) {}
Distribution::Distribution(std::string name, std::mt19937_64 &rng) : m_name(std::move(name)), m_rng(rng) {} BaseDistribution::BaseDistribution(std::string name, std::mt19937_64 &rng) : m_name(std::move(name)), m_rng(rng) {}
double Distribution::tau_wait() const { double BaseDistribution::tau_wait() const {
return std::exponential_distribution(1./m_tau_jump)(m_rng); return std::exponential_distribution(1./m_tau_jump)(m_rng);
} }
void Distribution::setParameters(const std::unordered_map<std::string, double> &parameters) { void BaseDistribution::setParameters(const std::unordered_map<std::string, double> &parameters) {
m_tau = parameters.at("tau"); m_tau = parameters.at("tau");
} }
std::unordered_map<std::string, double> Distribution::getParameters() const { std::unordered_map<std::string, double> BaseDistribution::getParameters() const {
return std::unordered_map<std::string, double>{ return std::unordered_map<std::string, double>{
{"tau", m_tau}, {"tau", m_tau},
}; };
} }
Distribution* Distribution::createFromInput(const std::string& input, std::mt19937_64& rng) { BaseDistribution* BaseDistribution::createFromInput(const std::string& input, std::mt19937_64& rng) {
if (input == "Delta") if (input == "Delta")
return new DeltaDistribution(rng); return new DeltaDistribution(rng);
@ -32,3 +33,4 @@ Distribution* Distribution::createFromInput(const std::string& input, std::mt199
throw std::invalid_argument("Invalid input " + input); throw std::invalid_argument("Invalid input " + input);
} }
}

View File

@ -4,12 +4,13 @@
#include <random> #include <random>
#include <unordered_map> #include <unordered_map>
class Distribution { namespace times {
class BaseDistribution {
public: public:
virtual ~Distribution() = default; virtual ~BaseDistribution() = default;
Distribution(std::string, double, std::mt19937_64&); BaseDistribution(std::string, double, std::mt19937_64&);
explicit Distribution(std::string, std::mt19937_64&); explicit BaseDistribution(std::string, std::mt19937_64&);
[[nodiscard]] double getTau() const { return m_tau; } [[nodiscard]] double getTau() const { return m_tau; }
void setTau(const double tau) { m_tau = tau; } void setTau(const double tau) { m_tau = tau; }
@ -24,7 +25,7 @@ public:
[[nodiscard]] virtual std::string toString() const = 0; [[nodiscard]] virtual std::string toString() const = 0;
static Distribution* createFromInput(const std::string& input, std::mt19937_64& rng); static BaseDistribution* createFromInput(const std::string& input, std::mt19937_64& rng);
protected: protected:
std::string m_name{"BaseDistribution"}; std::string m_name{"BaseDistribution"};
@ -32,5 +33,6 @@ protected:
double m_tau_jump{1.}; double m_tau_jump{1.};
std::mt19937_64& m_rng; std::mt19937_64& m_rng;
}; };
}
#endif //RWSIM_TIMESBASE_H #endif //RWSIM_TIMESBASE_H

View File

@ -1,7 +1,8 @@
#include "delta.h" #include "delta.h"
DeltaDistribution::DeltaDistribution(const double tau, std::mt19937_64& rng) : Distribution(std::string("Delta"), tau, rng) {} namespace times {
DeltaDistribution::DeltaDistribution(std::mt19937_64& rng) : Distribution(std::string("Delta"), rng) {} 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) {}
void DeltaDistribution::initialize() { void DeltaDistribution::initialize() {
m_tau_jump = m_tau; m_tau_jump = m_tau;
@ -12,4 +13,4 @@ void DeltaDistribution::draw_tau() {}
std::string DeltaDistribution::toString() const { std::string DeltaDistribution::toString() const {
return {"Delta/tau=" + std::to_string(m_tau)}; return {"Delta/tau=" + std::to_string(m_tau)};
} }
}

View File

@ -3,7 +3,8 @@
#include "base.h" #include "base.h"
class DeltaDistribution final : public Distribution { namespace times {
class DeltaDistribution final : public BaseDistribution {
public: public:
DeltaDistribution(double, std::mt19937_64&); DeltaDistribution(double, std::mt19937_64&);
explicit DeltaDistribution(std::mt19937_64 &rng); explicit DeltaDistribution(std::mt19937_64 &rng);
@ -13,5 +14,5 @@ public:
[[nodiscard]] std::string toString() const override; [[nodiscard]] std::string toString() const override;
}; };
}
#endif //RWSIM_TIMESDELTA_H #endif //RWSIM_TIMESDELTA_H

View File

@ -1,16 +1,17 @@
#include "lognormal.h" #include "lognormal.h"
#include <cmath> #include <cmath>
LogNormalDistribution::LogNormalDistribution(const double tau, const double sigma, std::mt19937_64& rng) : Distribution(std::string("LogNormal"), tau, rng), m_sigma(sigma), m_distribution(std::log(tau), sigma) {} namespace times {
LogNormalDistribution::LogNormalDistribution(std::mt19937_64& rng) : Distribution(std::string("LogNormal"), rng) {} 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) {}
void LogNormalDistribution::setParameters(const std::unordered_map<std::string, double> &parameters) { void LogNormalDistribution::setParameters(const std::unordered_map<std::string, double> &parameters) {
m_sigma = parameters.at("sigma"); m_sigma = parameters.at("sigma");
Distribution::setParameters(parameters); BaseDistribution::setParameters(parameters);
} }
std::unordered_map<std::string, double> LogNormalDistribution::getParameters() const { std::unordered_map<std::string, double> LogNormalDistribution::getParameters() const {
auto parameter = Distribution::getParameters(); auto parameter = BaseDistribution::getParameters();
parameter["sigma"] = m_sigma; parameter["sigma"] = m_sigma;
return parameter; return parameter;
} }
@ -27,3 +28,4 @@ void LogNormalDistribution::draw_tau() {
std::string LogNormalDistribution::toString() const { std::string LogNormalDistribution::toString() const {
return {"LogNormal/tau=" + std::to_string(m_tau) + "/sigma=" + std::to_string(m_sigma)}; return {"LogNormal/tau=" + std::to_string(m_tau) + "/sigma=" + std::to_string(m_sigma)};
} }
}

View File

@ -4,7 +4,8 @@
#include "base.h" #include "base.h"
#include <random> #include <random>
class LogNormalDistribution final : public Distribution { namespace times {
class LogNormalDistribution final : public BaseDistribution {
public: public:
LogNormalDistribution(double, double, std::mt19937_64&); LogNormalDistribution(double, double, std::mt19937_64&);
explicit LogNormalDistribution(std::mt19937_64 &rng); explicit LogNormalDistribution(std::mt19937_64 &rng);
@ -21,5 +22,6 @@ private:
double m_sigma{1}; double m_sigma{1};
std::lognormal_distribution<> m_distribution; std::lognormal_distribution<> m_distribution;
}; };
}
#endif //LOGNORMAL_H #endif //LOGNORMAL_H

View File

@ -1,9 +1,6 @@
#include <vector> #include <vector>
#include <array>
#include <chrono> #include <chrono>
#include <iostream> #include <iostream>
#include <cmath>
#include <utility>
int nearest_index(const std::vector<double> &x_ref, const double x, int start=0) { int nearest_index(const std::vector<double> &x_ref, const double x, int start=0) {

View File

@ -129,8 +129,8 @@ std::unordered_map<std::string, double> read_parameter(const std::filesystem::pa
std::string make_directory( std::string make_directory(
const Motion& motion, const motions::BaseMotion& motion,
const Distribution& distribution const times::BaseDistribution& distribution
) { ) {
std::ostringstream path_name; std::ostringstream path_name;
path_name << motion.toString() << "/" << distribution.toString(); path_name << motion.toString() << "/" << distribution.toString();

View File

@ -23,7 +23,7 @@ Arguments parse_args(int argc, char* argv[]);
std::pair<std::string, double> get_optional_parameter(std::vector<std::string>::const_iterator &it); std::pair<std::string, double> get_optional_parameter(std::vector<std::string>::const_iterator &it);
std::unordered_map<std::string, double> read_parameter(const std::filesystem::path&); std::unordered_map<std::string, double> read_parameter(const std::filesystem::path&);
std::string make_directory(const Motion&, const Distribution&); std::string make_directory(const motions::BaseMotion&, const times::BaseDistribution&);
void save_parameter_to_file(const std::string&, const std::string&, const std::unordered_map<std::string, double>&, const std::unordered_map<std::string, double>&); void save_parameter_to_file(const std::string&, const std::string&, const std::unordered_map<std::string, double>&, const std::unordered_map<std::string, double>&);
void save_data_to_file(const std::string&, const std::string&, const std::vector<double>&, const std::map<double, std::vector<double>>&, const std::unordered_map<std::string, double>&); void save_data_to_file(const std::string&, const std::string&, const std::vector<double>&, const std::map<double, std::vector<double>>&, const std::unordered_map<std::string, double>&);