NMRRWSims/cpp
1
0
Fork 0
Code Issues 3 Pull Requests Actions Packages Projects Releases Wiki Activity

Compare commits

base: NMRRWSims:master
Branches Tags
NMRRWSims:master
NMRRWSims:angle_dependence
..
compare: NMRRWSims:angle_dependence
Branches Tags
NMRRWSims:master
NMRRWSims:angle_dependence

2 Commits
master ... angle_depe

Author SHA1 Message Date
Dominik Demuth
78d008aa5a wip 2024-11-29 10:38:55 +01:00
Dominik Demuth
1bf927329d added flexibility 2024-11-28 19:40:07 +01:00
45 changed files with 553 additions and 1078 deletions
Show Stats Expand all files Collapse all files

61
angle_dependence.py
View File

@ -1,61 +0,0 @@
import pathlib
import re
import numpy as np
from matplotlib import pyplot as plt
from python.helpers import read_parameter_file
angles = []
tau_cc = []
tau_ss = []
tau2 = []
tevo = [5e-6, 10e-6]
for fit_files in pathlib.Path('.').glob(f'IsotropicAngle/angle=*/Delta/tau=*/ste_fit_*.dat'):
folder = fit_files.parent
file_base = fit_files.stem.replace('ste_fit_', '')
parameter = read_parameter_file(folder / ('ste_' + file_base + '_parameter.txt'))
angles.append(parameter['angle'])
with fit_files.open('r') as f:
# tau of F2 is hidden in the second header line
for _ in range(2):
line = f.readline()
tau2.append(float(re.search('tau=(.+?)\s', line).group(1)))
fit_values = np.loadtxt(fit_files)
x = fit_values[:, 0]
# get indexes for given evolution times
nearest_idx = [np.searchsorted(x, tt) for tt in tevo]
tau_cc.append(fit_values[nearest_idx, 1])
tau_ss.append(fit_values[nearest_idx, 4])
angles = np.array(angles)
tau_cc = np.array(tau_cc)
tau_ss = np.array(tau_ss)
tau2 = np.array(tau2)
sortidx = np.argsort(angles)
angles = angles[sortidx]
tau_cc = tau_cc[sortidx]
tau_ss = tau_ss[sortidx]
tau2 = tau2[sortidx]
plt.semilogy(angles, tau_cc, '-.')
plt.semilogy(angles, tau_ss, '--')
plt.semilogy(angles, tau2)
np.savetxt(
'tau_angles.dat',
np.c_[angles, tau_cc, tau_ss, tau2],
header=f"Angle dependence of correlation times\nfor evolution times {tevo}\nangle->tau_cc->tauss->tau2"
)
plt.show()

12
config.txt
View File

@ -1,21 +1,23 @@
# Simulation part # Simulation part
num_walker=20000 num_walker=20000
# Motion model part # Motion model part
delta=126e3 delta=161e3
eta=0.0 eta=0.0
# Distribution part # Distribution part
tau=1e-3 tau=1e-3
angle1=2
angle2=30
probability1=0
# Spectrum part # Spectrum part
dwell_time=1e-8 dwell_time=1e-8
num_acq=4096 num_acq=4096
techo_start=1e-6 techo_start=0e-6
techo_stop=40e-6 techo_stop=40e-6
techo_steps=5 techo_steps=5
# STE part # STE part
tevo_start=2e-6 tevo_start=2e-6
tevo_stop=120e-6 tevo_stop=50e-6
tevo_steps=121 tevo_steps=49
tmix_start=1e-5 tmix_start=1e-5
tmix_stop=1e1 tmix_stop=1e1
tmix_steps=31 tmix_steps=31
tpulse4=10e-6

80
main.py Normal file
View File

@ -0,0 +1,80 @@
import matplotlib.pyplot as plt
from python.ste import *
from python.helpers import *
# Simulation parameter
motion = 'IsotropicAngle'
distribution = 'Delta'
# parameter = {}
parameter = {
"angle": np.linspace(1, 180, num=180),
"eta": 0.1
# "sigma": 0.1,
# "tau": np.logspace(-1, 1, num=3)
}
parameter = prepare_rw_parameter(parameter)
# fig_tau_cc, ax_tau_cc = plt.subplots()
# ax_tau_cc.set_title('tau_cc')
#
# fig_beta_cc, ax_beta_cc = plt.subplots()
# ax_beta_cc.set_title('beta_cc')
#
# fig_finfty_cc, ax_finfty_cc = plt.subplots()
# ax_finfty_cc.set_title('f_infty_cc')
#
# fig_tau_ss, ax_tau_ss = plt.subplots()
# ax_tau_ss.set_title('tau_ss')
#
# fig_beta_ss, ax_beta_ss = plt.subplots()
# ax_beta_ss.set_title('beta_ss')
#
# fig_finfty_ss, ax_finfty_ss = plt.subplots()
# ax_finfty_ss.set_title('f_infty_ss')
tau_ss_angles = []
tau_cc_angles = []
tau_2_angles = []
angles = []
for variation in parameter:
print(f"\nRun RW for {motion}/{distribution} with arguments {variation}\n")
run_sims(motion, distribution, ste=True, spectrum=False, **variation)
conf_file = find_config_file(motion, distribution, variation)
vary_string, tau_cc, beta_cc, finfty_cc = fit_and_save_ste(conf_file, 'coscos', plot_decays=False, verbose=False)
_, tau_ss, beta_ss, finfty_ss = fit_and_save_ste(conf_file, 'sinsin', plot_decays=False, verbose=False)
_, tau_2, beta_2, finfty_2 = fit_and_save_ste(conf_file, 'f2', plot_decays=False, verbose=True)
# ax_tau_cc.semilogy(tau_cc[:, 0], tau_cc[:, 1], label=vary_string)
# ax_tau_cc.axhline(tau_2[:, 1], color='k', linestyle='--')
# ax_beta_cc.plot(*beta_cc.T, label=vary_string)
# ax_finfty_cc.plot(*finfty_cc.T, label=vary_string)
# ax_tau_ss.semilogy(tau_ss[:, 0], tau_ss[:, 1], label=vary_string)
# ax_tau_ss.axhline(tau_2[:, 1], color='k', linestyle='--')
# ax_beta_ss.plot(*beta_ss.T, label=vary_string)
# ax_finfty_ss.plot(*finfty_ss.T, label=vary_string)
angles.append(variation['angle'])
tau_ss_angles.append(tau_ss[0, 1])
tau_cc_angles.append(tau_cc[0, 1])
tau_2_angles.append(tau_2[0, 1])
fig, ax = plt.subplots()
ax.semilogy(angles, tau_ss_angles, 'o', label='SS (4.9mus)')
ax.plot(angles, tau_cc_angles, 'o', label='CC (4.9mus)')
ax.plot(angles, tau_2_angles, 'o', label='F2')
ax.legend()
np.savetxt('angle_eta.dat', np.c_[angles, tau_cc_angles, tau_ss_angles, tau_2_angles], header='#x\tcc\tss\tf2')
fig2, ax2 = plt.subplots()
ax2.plot(angles, np.array(tau_cc_angles)/np.array(tau_2_angles), 'o')
# for ax in [ax_tau_cc, ax_beta_cc, ax_finfty_cc, ax_tau_ss, ax_beta_ss, ax_finfty_ss]:
# ax.legend()
plt.show()

41
process_spectrum.py
View File

@ -1,41 +0,0 @@
import pathlib
import matplotlib.pyplot as plt
from python.spectrum import *
tpulse = 3e-6
gb = 4e3
fig_spec, ax_spec = plt.subplots()
# for conf_file in pathlib.Path('.').glob(f'SixSiteOctahedral/angle=59/Delta/tau=*/timesignal_*_parameter.txt'):
for conf_file in pathlib.Path('.').glob(f'RandomJumpOnCone/angle=128.*/Delta/tau=*/timesignal_*_parameter.txt'):
vary_string = post_process_spectrum(conf_file, tpulse=tpulse, apod=gb)
print(conf_file)
# ax_spec
#
# ax_tau_cc.semilogy(tau_cc[:, 0], tau_cc[:, 1], label=vary_string)
# ax_tau_cc.axhline(tau_2[:, 1], color='k', linestyle='--')
# ax_beta_cc.plot(*beta_cc.T, label=vary_string)
# ax_finfty_cc.plot(*finfty_cc.T, label=vary_string)
# ax_tau_ss.semilogy(tau_ss[:, 0], tau_ss[:, 1], label=vary_string)
# ax_tau_ss.axhline(tau_2[:, 1], color='k', linestyle='--')
# ax_beta_ss.plot(*beta_ss.T, label=vary_string)
# ax_finfty_ss.plot(*finfty_ss.T, label=vary_string)
#
# np.savetxt(
# conf_file.with_name(f'ste_fit_{vary_string}.dat'),
# np.c_[
# tau_cc, beta_cc[:, 1], finfty_cc[:, 1],
# tau_ss[:, 1], beta_ss[:, 1], finfty_ss[:, 1],
# ],
# header=f'Fit STE {vary_string}\n'
# f'F2: tau={tau_2[0, 1]} beta={beta_2[0, 1]} finfty={finfty_2[0, 1]}\n'
# f'tevo\ttaucc\tbetacc\tfinftycc\ttauss\tbetass\tfinftyss',
# )
#
# for ax in [ax_tau_cc, ax_beta_cc, ax_finfty_cc, ax_tau_ss, ax_beta_ss, ax_finfty_ss]:
# ax.legend()
# plt.show()

49
process_ste.py
View File

@ -1,49 +0,0 @@
from python.ste import *
fig_tau_cc, ax_tau_cc = plt.subplots()
ax_tau_cc.set_title('tau_cc')
fig_beta_cc, ax_beta_cc = plt.subplots()
ax_beta_cc.set_title('beta_cc')
fig_finfty_cc, ax_finfty_cc = plt.subplots()
ax_finfty_cc.set_title('f_infty_cc')
fig_tau_ss, ax_tau_ss = plt.subplots()
ax_tau_ss.set_title('tau_ss')
fig_beta_ss, ax_beta_ss = plt.subplots()
ax_beta_ss.set_title('beta_ss')
fig_finfty_ss, ax_finfty_ss = plt.subplots()
ax_finfty_ss.set_title('f_infty_ss')
for conf_file in pathlib.Path('.').glob(f'IsotropicAngle/angle=*/Delta/tau=*/*_parameter.txt'):
print(conf_file)
vary_string, tau_cc, beta_cc, finfty_cc = fit_ste(conf_file, f'coscos', plot_decays=False, verbose=False)
_, tau_ss, beta_ss, finfty_ss = fit_ste(conf_file, f'sinsin', plot_decays=False, verbose=False)
_, tau_2, beta_2, finfty_2 = fit_ste(conf_file, f'f2', plot_decays=True, verbose=True)
ax_tau_cc.semilogy(tau_cc[:, 0], tau_cc[:, 1], label=vary_string)
ax_tau_cc.axhline(tau_2[:, 1], color='k', linestyle='--')
ax_beta_cc.plot(*beta_cc.T, label=vary_string)
ax_finfty_cc.plot(*finfty_cc.T, label=vary_string)
ax_tau_ss.semilogy(tau_ss[:, 0], tau_ss[:, 1], label=vary_string)
ax_tau_ss.axhline(tau_2[:, 1], color='k', linestyle='--')
ax_beta_ss.plot(*beta_ss.T, label=vary_string)
ax_finfty_ss.plot(*finfty_ss.T, label=vary_string)
np.savetxt(
conf_file.with_name(f'ste_fit_{vary_string}.dat'),
np.c_[
tau_cc, beta_cc[:, 1], finfty_cc[:, 1],
tau_ss[:, 1], beta_ss[:, 1], finfty_ss[:, 1],
],
header=f'Fit STE {vary_string}\n'
f'F2: tau={tau_2[0, 1]} beta={beta_2[0, 1]} finfty={finfty_2[0, 1]}\n'
f'tevo\ttaucc\tbetacc\tfinftycc\ttauss\tbetass\tfinftyss',
)
for ax in [ax_tau_cc, ax_beta_cc, ax_finfty_cc, ax_tau_ss, ax_beta_ss, ax_finfty_ss]:
ax.legend()
plt.show()

25
python/helpers.py
View File

@ -1,10 +1,10 @@
from __future__ import annotations from __future__ import annotations
import pathlib import pathlib
import re
import subprocess import subprocess
from itertools import product from itertools import product
def prepare_rw_parameter(parameter: dict) -> list: def prepare_rw_parameter(parameter: dict) -> list:
""" """
Converts a dictionary of iterables to list of dictionaries Converts a dictionary of iterables to list of dictionaries
@ -53,11 +53,23 @@ def run_sims(
subprocess.run(arguments) subprocess.run(arguments)
def find_config_file(config_path: str | pathlib.Path, varied_params: dict[str, float]) -> dict[str, float]: def find_config_file(motion: str, distribution: str, var_params: dict) -> pathlib.Path:
parameter = read_parameter_file(config_path) # TODO handle situation if multiple files fit
parameter.update(varied_params) p_file = None
if var_params:
var_string = '|'.join(([f'{k}={v:1.6e}' for (k, v) in var_params.items()])).replace('.', '\.').replace('+', '\+')
pattern = re.compile(var_string)
for p_file in pathlib.Path('.').glob('*_parameter.txt'):
if len(re.findall(pattern, str(p_file))) == len(var_params) and re.search(f'{motion}_{distribution}', str(p_file)):
return p_file
raise ValueError(f'No parameter file found for {motion}, {distribution}, {var_params}')
else:
for p_file in pathlib.Path('.').glob('*_parameter.txt'):
if re.search(f'{motion}_{distribution}', str(p_file)):
return p_file
raise ValueError(f'No parameter file found for {motion}, {distribution}, {var_params}')
return parameter
def read_parameter_file(path: str | pathlib.Path) -> dict[str, float]: def read_parameter_file(path: str | pathlib.Path) -> dict[str, float]:
@ -68,9 +80,8 @@ def read_parameter_file(path: str | pathlib.Path) -> dict[str, float]:
parameter_dict = {} parameter_dict = {}
with path.open('r') as f: with path.open('r') as f:
for line in f.readlines(): for line in f.readlines():
if line.startswith('#'):
continue
k, v = line.split('=') k, v = line.split('=')
parameter_dict[k] = float(v) parameter_dict[k] = float(v)
k, v = line.split('=')
return parameter_dict return parameter_dict

52
python/spectrum.py
View File

@ -1,7 +1,8 @@
import numpy
import numpy as np import numpy as np
import matplotlib.pyplot as plt from matplotlib import pyplot
from python.helpers import read_parameter_file
# parameter for spectrum simulations # parameter for spectrum simulations
lb = 2e3 lb = 2e3
@ -40,39 +41,30 @@ def pulse_attn(freq: np.ndarray, t_pulse: float) -> np.ndarray:
return np.pi * numerator / denominator / 2 return np.pi * numerator / denominator / 2
def post_process_spectrum(parameter_file, apod, tpulse): def post_process_spectrum(taus, apod, tpulse):
parameter = read_parameter_file(parameter_file) reduction_factor = np.zeros((taus.size, 5)) # hard-coded t_echo :(
# files have form ste_arg=0.000000e+01_parameter.txt, first remove ste part then parameter.txt to get variables for i, tau in enumerate(taus):
varied_string = str(parameter_file).partition('_')[-1].rpartition('_')[0] try:
raw_data = np.loadtxt(f'fid_tau={tau:.6e}.dat')
except OSError:
continue
# make evolution times t = raw_data[:, 0]
tevo = np.linspace(parameter['techo_start'], parameter['techo_stop'], num=int(parameter['techo_steps'])) timesignal = raw_data[:, 1:]
if varied_string: timesignal *= dampening(t, apod)[:, None]
raw_data = np.loadtxt(parameter_file.with_name(f'timesignal_{varied_string}.dat')) timesignal[0, :] /= 2
else:
raw_data = np.loadtxt(parameter_file.with_name(f'timesignal.dat'))
t = raw_data[:, 0] # FT to spectrum
timesignal = raw_data[:, 1:] freq = np.fft.fftshift(np.fft.fftfreq(t.size, d=1e-6))
spec = np.fft.fftshift(np.fft.fft(timesignal, axis=0), axes=0).real
spec *= pulse_attn(freq, t_pulse=tpulse)[:, None]
timesignal *= dampening(t, apod)[:, None] reduction_factor[i, :] = 2*timesignal[0, :]
timesignal[0, :] /= 2
# FT to spectrum plt.plot(freq, spec)
freq = np.fft.fftshift(np.fft.fftfreq(t.size, d=parameter['dwell_time'])) plt.show()
spec = np.fft.fftshift(np.fft.fft(timesignal, axis=0), axes=0).real
spec *= pulse_attn(freq, t_pulse=tpulse)[:, None]
# plt.semilogx(taus, reduction_factor, '.')
#
# reduction_factor[i, :] = 2*timesignal[0, :]
plt.plot(freq, spec)
plt.gca().set_title(varied_string)
plt.show() plt.show()
#
# plt.semilogx(taus, reduction_factor, '.')
# plt.show()

15
python/ste.py
View File

@ -1,5 +1,6 @@
import pathlib import pathlib
import numpy
import numpy as np import numpy as np
from matplotlib import pyplot as plt from matplotlib import pyplot as plt
from scipy.optimize import curve_fit from scipy.optimize import curve_fit
@ -58,12 +59,12 @@ def fit_decay(x: np.ndarray, y: np.ndarray, tevo: np.ndarray, verbose: bool = Tr
return tau_fit, beta_fit, finfty_fit return tau_fit, beta_fit, finfty_fit
def fit_ste( def fit_and_save_ste(
parameter_file: pathlib.Path, parameter_file: pathlib.Path,
prefix: str, prefix: str,
plot_decays: bool = True, plot_decays: bool = True,
verbose: bool = True, verbose: bool = True,
) -> tuple[str, np.ndarray, np.ndarray, np.ndarray]: ) -> tuple[str, np.ndarray, np.ndarray, np.ndarray]:
# read simulation parameters # read simulation parameters
parameter = read_parameter_file(parameter_file) parameter = read_parameter_file(parameter_file)
@ -73,10 +74,7 @@ def fit_ste(
# make evolution times # make evolution times
tevo = np.linspace(parameter['tevo_start'], parameter['tevo_stop'], num=int(parameter['tevo_steps'])) tevo = np.linspace(parameter['tevo_start'], parameter['tevo_stop'], num=int(parameter['tevo_steps']))
if varied_string: raw_data = np.loadtxt(f'{prefix}_{varied_string}.dat')
raw_data = np.loadtxt(parameter_file.with_name(f'{prefix}_{varied_string}.dat'))
else:
raw_data = np.loadtxt(parameter_file.with_name(f'{prefix}.dat'))
t_mix = raw_data[:, 0] t_mix = raw_data[:, 0]
decay = raw_data[:, 1:] decay = raw_data[:, 1:]
@ -91,4 +89,9 @@ def fit_ste(
print(f'Fit {prefix}') print(f'Fit {prefix}')
tau, beta, finfty = fit_decay(t_mix, decay, tevo, verbose=verbose) tau, beta, finfty = fit_decay(t_mix, decay, tevo, verbose=verbose)
np.savetxt(f'tau_{prefix}_{varied_string}.dat', tau)
np.savetxt(f'beta_{prefix}_{varied_string}.dat', beta)
np.savetxt(f'finfty_{prefix}_{varied_string}.dat', finfty)
return varied_string, tau, beta, finfty return varied_string, tau, beta, finfty

39
src/CMakeLists.txt
View File

@ -1,16 +1,35 @@
cmake_minimum_required(VERSION 3.18)
add_subdirectory(times) set(CMAKE_CXX_STANDARD 17)
add_subdirectory(motions)
add_subdirectory(utils)
add_library(simulation STATIC sims.cpp sims.h)
target_link_libraries(simulation PRIVATE utils)
add_executable( add_executable(rwsim main.cpp
rwsim utils/functions.h
main.cpp utils/functions.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)

7
src/main.cpp
View File

@ -1,9 +1,10 @@
#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>
@ -36,8 +37,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());
motions::BaseMotion *motion = motions::BaseMotion::createFromInput(args.motion_type, rng); Motion *motion = Motion::createFromInput(args.motion_type, rng);
times::BaseDistribution *dist = times::BaseDistribution::createFromInput(args.distribution_type, rng); Distribution *dist = Distribution::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);
} }

25
src/motions/CMakeLists.txt
View File

@ -1,25 +0,0 @@
# 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
nsiteconejump.cpp
nsiteconejump.h
)

133
src/motions/base.cpp
View File

@ -1,89 +1,64 @@
#include "base.h" #include "base.h"
#include <iostream>
#include "coordinates.h" #include "coordinates.h"
#include "bimodalangle.h" #include "bimodalangle.h"
#include "isosmallangle.h" #include "isosmallangle.h"
#include "random.h" #include "random.h"
#include "foursitejump.h" #include "tetrahedral.h"
#include "nsiteconejump.h"
#include "sixsitejump.h"
#include "rjoac.h"
#include <stdexcept> #include <stdexcept>
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) {
m_uni_dist = std::uniform_real_distribution(0., 1.);
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, std::mt19937_64& rng) : m_name(std::move(name)), m_rng(rng) {
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;
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));
}
double BaseMotion::omega_q(const coordinates::SphericalPos& pos) const {
auto [cos_theta, phi] = pos;
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);
return {cos_theta, phi};
}
BaseMotion* BaseMotion::createFromInput(const std::string& input, std::mt19937_64& rng) {
if (input == "FourSiteTetrahedral")
return new FourSiteTetrahedron(rng);
if (input == "SixSiteOctahedralC3")
return new SixSiteOctahedronC3(rng);
if (input == "IsotropicAngle")
return new SmallAngle(rng);
if (input == "RandomJump")
return new RandomJump(rng);
if (input == "BimodalAngle")
return new BimodalAngle(rng);
if (input == "NSiteConeJump")
return new NSiteJumpOnCone(rng);
if (input == "RandomJumpOnCone")
return new RandomJumpOnCone(rng);
throw std::invalid_argument("Invalid input " + input);
}
void BaseMotion::setParameters(const std::unordered_map<std::string, double> &parameters) {
m_delta = parameters.at("delta");
m_eta = parameters.at("eta");
}
std::unordered_map<std::string, double> BaseMotion::getParameters() const {
return std::unordered_map<std::string, double>{
{"delta", m_delta},
{"eta", m_eta}
};
}
std::ostream& operator<<(std::ostream& os, const BaseMotion& m) {
os << m.getName();
return os;
}
} }
Motion::Motion(std::string name, std::mt19937_64& rng) : m_name(std::move(name)), m_rng(rng) {
m_uni_dist = std::uniform_real_distribution(0., 1.);
}
double Motion::omega_q(const double cos_theta, const double phi) const {
const double cos_theta_square = cos_theta * cos_theta;
const double sin_theta_square = 1. - cos_theta_square;
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 {
auto [cos_theta, phi] = pos;
return omega_q(cos_theta, phi);
}
SphericalPos Motion::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);
return {cos_theta, phi};
}
Motion* Motion::createFromInput(const std::string& input, std::mt19937_64& rng) {
if (input == "TetrahedralJump")
return new TetrahedralJump(rng);
if (input == "IsotropicAngle")
return new SmallAngle(rng);
if (input == "RandomJump")
return new RandomJump(rng);
if (input == "BimodalAngle")
return new BimodalAngle(rng);
throw std::invalid_argument("Invalid input " + input);
}
void Motion::setParameters(const std::unordered_map<std::string, double> &parameters) {
m_delta = parameters.at("delta");
m_eta = parameters.at("eta");
}
std::ostream& operator<<(std::ostream& os, const Motion& m) {
os << m.getName();
return os;
}

59
src/motions/base.h
View File

@ -6,43 +6,40 @@
#include <random> #include <random>
#include <unordered_map> #include <unordered_map>
namespace motions { class Motion {
class BaseMotion { public:
public: virtual ~Motion() = default;
virtual ~BaseMotion() = default;
BaseMotion(std::string, double, double, std::mt19937_64&); Motion(std::string, double, double, std::mt19937_64&);
explicit BaseMotion(std::string, std::mt19937_64&); explicit Motion(std::string, std::mt19937_64&);
coordinates::SphericalPos draw_position(); SphericalPos draw_position();
[[nodiscard]] double omega_q(double, double) const; [[nodiscard]] double omega_q(double, double) const;
[[nodiscard]] double omega_q(const coordinates::SphericalPos&) const; [[nodiscard]] double omega_q(const SphericalPos&) const;
virtual void initialize() = 0; virtual void initialize() = 0;
virtual double jump() = 0; virtual double jump() = 0;
virtual void setParameters(const std::unordered_map<std::string, double>&); virtual void setParameters(const std::unordered_map<std::string, double>&);
[[nodiscard]] virtual std::unordered_map<std::string, double> getParameters() const;
[[nodiscard]] double getDelta() const { return m_delta; }
void setDelta(const double delta) { m_delta = delta; }
[[nodiscard]] double getEta() const { return m_eta; }
void setEta(const double eta) { m_eta = eta; }
[[nodiscard]] std::string getName() const { return m_name; }
[[nodiscard]] double getInitOmega() const { return m_initial_omega; }
[[nodiscard]] virtual std::string toString() const = 0; [[nodiscard]] double getDelta() const { return m_delta; }
void setDelta(const double delta) { m_delta = delta; }
[[nodiscard]] double getEta() const { return m_eta; }
void setEta(const double eta) { m_eta = eta; }
[[nodiscard]] std::string getName() const { return m_name; }
[[nodiscard]] double getInitOmega() const { return m_initial_omega; };
static BaseMotion* createFromInput(const std::string& input, std::mt19937_64& rng); static Motion* createFromInput(const std::string& input, std::mt19937_64& rng);
protected: protected:
std::string m_name{"BaseMotion"}; std::string m_name{"BaseMotion"};
double m_delta{1.}; double m_delta{1.};
double m_eta{0.}; double m_eta{0.};
std::mt19937_64& m_rng; std::mt19937_64& m_rng;
std::uniform_real_distribution<> m_uni_dist; std::uniform_real_distribution<> m_uni_dist;
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

61
src/motions/bimodalangle.cpp
View File

@ -1,46 +1,31 @@
#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();
m_initial_omega = omega_q(m_prev_pos); m_initial_omega = omega_q(m_prev_pos);
} };
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 = 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) {
BaseMotion::setParameters(parameter); Motion::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.;
m_prob = parameter.at("probability1"); m_prob = parameter.at("probability1");
}
std::unordered_map<std::string, double> BimodalAngle::getParameters() const {
auto parameter = BaseMotion::getParameters();
parameter["angle1"] = m_angle1 * 180 / M_PI;
parameter["angle2"] = m_angle2 * 180 / M_PI;
parameter["probality1"] = m_prob;
return parameter;
}
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)};
}
} }

32
src/motions/bimodalangle.h
View File

@ -3,25 +3,21 @@
#define BIMODALANGLE_H #define BIMODALANGLE_H
#include "base.h" #include "base.h"
#include "coordinates.h"
namespace motions { class BimodalAngle : public Motion {
class BimodalAngle final : 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&);
void initialize() override; void initialize() override;
double jump() override; double jump() override;
void setParameters(const std::unordered_map<std::string, double> &) 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; protected:
double m_angle1{0};
double m_angle2{0};
double m_prob{0};
SphericalPos m_prev_pos{0., 0.};
};
protected:
double m_angle1{0};
double m_angle2{0};
double m_prob{0};
coordinates::SphericalPos m_prev_pos{0., 0.};
};
}
#endif //BIMODALANGLE_H #endif //BIMODALANGLE_H

37
src/motions/conewobble.cpp
View File

@ -1,37 +0,0 @@
#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(std::mt19937_64 &rng) : BaseMotion("Wobble in Cone", rng) {}
void WobbleCone::initialize() {
m_axis = draw_position();
}
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);
return omega_q(rotate(m_axis, real_angle, phi));
}
void WobbleCone::setParameters(const std::unordered_map<std::string, double> &parameters) {
BaseMotion::setParameters(parameters);
m_angle = parameters.at("angle");
}
std::unordered_map<std::string, double> WobbleCone::getParameters() const {
auto parameter = BaseMotion::getParameters();
parameter["angle"] = m_angle;
return parameter;
}
std::string WobbleCone::toString() const {
return std::string("ConeWobble/angle=") + std::to_string(m_angle);
}
}

30
src/motions/conewobble.h
View File

@ -1,30 +0,0 @@
#ifndef CONEWOBBLE_H
#define CONEWOBBLE_H
#include "base.h"
#include "coordinates.h"
#include <random>
namespace motions {
class WobbleCone final: public BaseMotion {
public:
WobbleCone(double, double, double, std::mt19937_64&);
explicit WobbleCone(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 //CONEWOBBLE_H

58
src/motions/coordinates.cpp
View File

@ -3,39 +3,37 @@
#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)};
const double sin_beta{std::sin(beta)}; const double sin_beta{std::sin(beta)};
const double cos_beta{std::cos(beta)}; const double cos_beta{std::cos(beta)};
const double cos_theta{old_pos.cos_theta}; const double cos_theta{old_pos.cos_theta};
if (std::abs(cos_theta) == 1) { if (std::abs(cos_theta) == 1) {
return xyz_to_spherical(CartesianPos{cos_alpha * cos_beta, cos_alpha * sin_beta, cos_alpha * cos_theta}); return xyz_to_spherical(CartesianPos{cos_alpha * cos_beta, cos_alpha * sin_beta, cos_alpha * cos_theta});
}
const double norm{std::sqrt(1 - cos_theta * cos_theta) + 1e-15};
auto [x, y , z] = spherical_to_xyz(old_pos);
const auto new_pos = CartesianPos{
cos_alpha * x + sin_alpha * (-cos_beta * y - sin_beta * x * z) / norm,
cos_alpha * y + sin_alpha * (cos_beta * x - sin_beta * y * z) / norm,
cos_alpha * z + sin_alpha * sin_beta * norm,
};
return xyz_to_spherical(new_pos);
} }
CartesianPos spherical_to_xyz(const SphericalPos& pos) { const double norm{std::sqrt(1 - cos_theta * cos_theta) + 1e-15};
const double sin_theta = std::sin(std::acos(pos.cos_theta));
return {sin_theta * std::cos(pos.phi), sin_theta * std::sin(pos.phi), pos.cos_theta};
}
SphericalPos xyz_to_spherical(const CartesianPos& pos) { auto [x, y , z] = spherical_to_xyz(old_pos);
return {pos.z, std::atan2(pos.y, pos.x)};
} const auto new_pos = CartesianPos{
} cos_alpha * x + sin_alpha * (-x * z * sin_beta - y * cos_beta) / norm,
cos_alpha * y + sin_alpha * (-y * z * sin_beta + x * cos_beta) / norm,
cos_alpha * z + sin_alpha * norm * sin_beta
};
return xyz_to_spherical(new_pos);
}
CartesianPos spherical_to_xyz(const SphericalPos& pos) {
const double sin_theta = std::sin(std::acos(pos.cos_theta));
return {sin_theta * std::cos(pos.phi), sin_theta * std::sin(pos.phi), pos.cos_theta};
}
SphericalPos xyz_to_spherical(const CartesianPos& pos) {
return {pos.z, std::atan2(pos.y, pos.x)};
}

26
src/motions/coordinates.h
View File

@ -2,21 +2,19 @@
#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; };
};
struct CartesianPos { struct CartesianPos {
double x; double x;
double y; double y;
double z; double z;
}; };
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

35
src/motions/foursitejump.cpp
View File

@ -1,35 +0,0 @@
#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 = 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"};
}
}

31
src/motions/foursitejump.h
View File

@ -1,31 +0,0 @@
#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

43
src/motions/isosmallangle.cpp
View File

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

27
src/motions/isosmallangle.h
View File

@ -5,21 +5,18 @@
#include "base.h" #include "base.h"
#include "coordinates.h" #include "coordinates.h"
namespace motions { class SmallAngle final : public Motion {
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&);
void initialize() override; void initialize() override;
double jump() override; double jump() override;
void setParameters(const std::unordered_map<std::string, double> &) 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_chi{0};
SphericalPos m_prev_pos{0., 0.};
};
private:
double m_chi{0};
coordinates::SphericalPos m_prev_pos{0., 0.};
};
}
#endif //RWSIM_MOTIONISOSMALLANGLE_H #endif //RWSIM_MOTIONISOSMALLANGLE_H

59
src/motions/nsiteconejump.cpp
View File

@ -1,59 +0,0 @@
//
// 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) :
BaseMotion("NSiteJumpOnCone", delta, eta, rng),
m_num_sites(num_sites),
m_chi(chi*M_PI/180.) {}
NSiteJumpOnCone::NSiteJumpOnCone(std::mt19937_64 &rng) : BaseMotion("NSiteJumpOnCone", rng) { }
void NSiteJumpOnCone::initialize() {
m_sites = std::vector<double>(m_num_sites);
m_chooser = std::uniform_int_distribution<>{1, m_num_sites - 1};
m_axis = draw_position();
const double alpha = m_uni_dist(m_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);
m_cone_idx %= m_num_sites;
return m_sites[m_cone_idx];
}
void NSiteJumpOnCone::setParameters(const std::unordered_map<std::string, double> &parameters) {
BaseMotion::setParameters(parameters);
m_num_sites = static_cast<int>(parameters.at("num_sites"));
m_chi = parameters.at("chi") * M_PI/180.;
}
std::unordered_map<std::string, double> NSiteJumpOnCone::getParameters() const {
auto parameter = BaseMotion::getParameters();
parameter["num_sites"] = m_num_sites;
parameter["chi"] = m_chi * 180. / M_PI;
return parameter;
}
std::string NSiteJumpOnCone::toString() const {
return std::to_string(m_num_sites) + "SiteJumpOnCone/angle=" + std::to_string(m_chi*180/M_PI);
}
} // motions

33
src/motions/nsiteconejump.h
View File

@ -1,33 +0,0 @@
#ifndef NSITEJUMPONCONE_H
#define NSITEJUMPONCONE_H
#include "base.h"
#include <random>
#include <vector>
namespace motions {
class NSiteJumpOnCone final : public BaseMotion {
public:
NSiteJumpOnCone(double, double, int, double, std::mt19937_64&);
explicit NSiteJumpOnCone(std::mt19937_64&);
void initialize() override;
double jump() override;
[[nodiscard]] std::string toString() const override;
void setParameters(const std::unordered_map<std::string, double> &) override;
[[nodiscard]] std::unordered_map<std::string, double> getParameters() const override;
private:
int m_num_sites{2};
std::vector<double> m_sites{};
double m_chi{1./2.};
int m_cone_idx = 0;
coordinates::SphericalPos m_axis{1, 0};
std::uniform_int_distribution<> m_chooser;
};
} // motions
#endif //NSITEJUMPONCONE_H

21
src/motions/random.cpp
View File

@ -1,20 +1,15 @@
#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(std::mt19937_64 &rng) : BaseMotion(std::string("RandomJump"), rng) {} RandomJump::RandomJump(const double delta, const double eta, std::mt19937_64 &rng) : Motion(std::string("RandomJump"), delta, eta, rng) {}
std::string RandomJump::toString() const { RandomJump::RandomJump(std::mt19937_64 &rng) : Motion(std::string("RandomJump"), rng) {}
return {"RandomJump"};
}
void RandomJump::initialize() { void RandomJump::initialize() {
m_initial_omega = RandomJump::jump(); m_initial_omega = RandomJump::jump();
} }
double RandomJump::jump() { double RandomJump::jump() {
return omega_q(draw_position()); return omega_q(draw_position());
}
} }

18
src/motions/random.h
View File

@ -5,17 +5,13 @@
#include "base.h" #include "base.h"
#include <random> #include <random>
namespace motions { class RandomJump final : public Motion {
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&);
[[nodiscard]] std::string toString() const override; void initialize() override;
double jump() override;
void initialize() override; };
double jump() override;
};
}
#endif //RWSIM_MOTIONRANDOMJUMP_H #endif //RWSIM_MOTIONRANDOMJUMP_H

32
src/motions/rjoac.cpp
View File

@ -1,32 +0,0 @@
#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(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
View File

@ -1,31 +0,0 @@
#ifndef RJOAC_H
#define RJOAC_H
#include "base.h"
#include "coordinates.h"
#include <random>
namespace motions {
class RandomJumpOnCone final: 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

41
src/motions/sixsitejump.cpp
View File

@ -1,41 +0,0 @@
#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) :
m_chi{chi*M_PI/180.},
BaseMotion(std::string{"SixSiteOctahedral"}, delta, eta, rng) {}
SixSiteOctahedronC3::SixSiteOctahedronC3(std::mt19937_64& rng) : BaseMotion(std::string{"SixSiteOctahedralC3"}, rng) {}
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);
const double m_chi_opposite = M_PI - m_chi;
for (int i = 0; i<3; i++) {
m_corners[2*i] = omega_q(rotate(c3_axis, m_chi, alpha + i * 2./3.*M_PI));
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();
}
double SixSiteOctahedronC3::jump() {
m_corner_idx += m_chooser(m_rng);
m_corner_idx %= 6;
return m_corners[m_corner_idx];
}
std::string SixSiteOctahedronC3::toString() const {
return {"SixSiteOctahedral/angle=" + std::to_string(m_chi / M_PI * 180.)};
}
}

32
src/motions/sixsitejump.h
View File

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

30
src/motions/tetrahedral.cpp Normal file
View File

@ -0,0 +1,30 @@
#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];
}

28
src/motions/tetrahedral.h Normal file
View File

@ -0,0 +1,28 @@
#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;
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

52
src/sims.cpp → src/simulation/sims.cpp
View File

@ -1,8 +1,9 @@
#include "sims.h" #include "sims.h"
#include "times/base.h" #include "../motions/base.h"
#include "utils/functions.h" #include "../times/base.h"
#include "utils/ranges.h" #include "../utils/functions.h"
#include "utils/io.h" #include "../utils/ranges.h"
#include "../utils/io.h"
#include <iostream> #include <iostream>
#include <algorithm> #include <algorithm>
@ -14,12 +15,11 @@
#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,
motions::BaseMotion& motion, Motion& motion,
times::BaseDistribution& dist Distribution& dist
) { ) {
const int num_walker = static_cast<int>(parameter["num_walker"]); const int num_walker = static_cast<int>(parameter["num_walker"]);
@ -72,9 +72,8 @@ void run_spectrum(
} }
// write fid to files // write fid to files
const auto path = make_directory(motion, dist); save_parameter_to_file("timesignal", motion.getName(), dist.getName(), parameter, optional);
save_parameter_to_file(std::string("timesignal"), path, parameter, optional); save_data_to_file("timesignal", motion.getName(), dist.getName(), t_fid, fid_dict, optional);
save_data_to_file(std::string("timesignal"), path, t_fid, fid_dict, optional);
printEnd(start); printEnd(start);
} }
@ -83,15 +82,14 @@ 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,
motions::BaseMotion& motion, Motion& motion,
times::BaseDistribution& dist Distribution& 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")]);
const int num_mix_times = static_cast<int>(parameter["tmix_steps"]); const int num_mix_times = static_cast<int>(parameter[std::string("tmix_steps")]);
const std::vector<double> evolution_times = linspace(parameter["tevo_start"], parameter["tevo_stop"], static_cast<int>(parameter["tevo_steps"])); const std::vector<double> evolution_times = linspace(parameter["tevo_start"], parameter["tevo_stop"], static_cast<int>(parameter["tevo_steps"]));
const std::vector<double> mixing_times = logspace(parameter["tmix_start"], parameter["tmix_stop"], num_mix_times); const std::vector<double> mixing_times = logspace(parameter["tmix_start"], parameter["tmix_stop"], num_mix_times);
const double tpulse4 = parameter["tpulse4"];
// make ste decay vectors and set them to zero // make ste decay vectors and set them to zero
std::map<double, std::vector<double>> cc_dict; std::map<double, std::vector<double>> cc_dict;
@ -105,7 +103,7 @@ void run_ste(
std::vector<double> f2(num_mix_times); std::vector<double> f2(num_mix_times);
// each trajectory must have a duration of at least tmax // each trajectory must have a duration of at least tmax
const double tmax = *std::max_element(evolution_times.begin(), evolution_times.end()) * 2 + *std::max_element(mixing_times.begin(), mixing_times.end()) + 2*tpulse4; const double tmax = *std::max_element(evolution_times.begin(), evolution_times.end()) * 2 + *std::max_element(mixing_times.begin(), mixing_times.end());
// set parameter in distribution and motion model // set parameter in distribution and motion model
dist.setParameters(parameter); dist.setParameters(parameter);
@ -129,7 +127,6 @@ void run_ste(
f2[f2_idx] += traj_omega[f2_pos] * motion.getInitOmega() / num_walker; f2[f2_idx] += traj_omega[f2_pos] * motion.getInitOmega() / num_walker;
} }
for (auto& [t_evo_j, _] : cc_dict) { for (auto& [t_evo_j, _] : cc_dict) {
auto& cc_j = cc_dict[t_evo_j]; auto& cc_j = cc_dict[t_evo_j];
auto& ss_j = ss_dict[t_evo_j]; auto& ss_j = ss_dict[t_evo_j];
@ -141,20 +138,16 @@ void run_ste(
const double ss_tevo = std::sin(dephased); const double ss_tevo = std::sin(dephased);
for (int mix_idx = 0; mix_idx < num_mix_times; mix_idx++) { for (int mix_idx = 0; mix_idx < num_mix_times; mix_idx++) {
// get phase at end of mixing time // get phase at end of mixing time
const double time_end_mix = mixing_times[mix_idx] + t_evo_j; 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);
const double phase_4pulse = lerp(traj_time, traj_phase, time_pulse4, current_pos);
// get phase at echo position // get phase at echo position
const double time_echo = time_pulse4 + tpulse4 + t_evo_j; const double time_echo = mixing_times[mix_idx] + 2 * 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; const double rephased = lerp(traj_time, traj_phase, time_echo, current_pos) - phase_mix_end;
cc_j[mix_idx] += cc_tevo * std::cos(rephased) / num_walker; cc_j[mix_idx] += cc_tevo * std::cos(rephased) / num_walker;
ss_j[mix_idx] += ss_tevo * std::sin(rephased) / num_walker; ss_j[mix_idx] += ss_tevo * std::sin(rephased) / num_walker;
@ -164,19 +157,18 @@ void run_ste(
} }
// write to files // write to files
const auto folders = make_directory(motion, dist); save_parameter_to_file("ste", motion.getName(), dist.getName(), parameter, optional);
save_parameter_to_file(std::string("ste"), folders, parameter, optional); save_data_to_file("coscos", motion.getName(), dist.getName(), mixing_times, cc_dict, optional);
save_data_to_file(std::string("coscos"), folders, mixing_times, cc_dict, optional); save_data_to_file("sinsin", motion.getName(), dist.getName(), mixing_times, ss_dict, optional);
save_data_to_file(std::string("sinsin"), folders, mixing_times, ss_dict, optional); save_data_to_file("f2", motion.getName(), dist.getName(), mixing_times, f2, optional);
save_data_to_file(std::string("f2"), folders, mixing_times, f2, optional);
printEnd(start); printEnd(start);
} }
void make_trajectory( void make_trajectory(
motions::BaseMotion& motion, Motion& motion,
times::BaseDistribution& dist, Distribution& 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,

10
src/sims.h → src/simulation/sims.h
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, motions::BaseMotion& motion, times::BaseDistribution& dist); void run_spectrum(std::unordered_map<std::string, double>& parameter, std::unordered_map<std::string, double> optional, Motion& motion, Distribution& 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, motions::BaseMotion& motion, times::BaseDistribution& dist); void run_ste(std::unordered_map<std::string, double>& parameter, std::unordered_map<std::string, double> optional, Motion& motion, Distribution& 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(motions::BaseMotion& motion, times::BaseDistribution& dist, double t_max, std::vector<double>& out_time, std::vector<double>& out_phase, std::vector<double>& out_omega); 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);
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);

37
src/times/base.cpp
View File

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

43
src/times/base.h
View File

@ -4,35 +4,30 @@
#include <random> #include <random>
#include <unordered_map> #include <unordered_map>
namespace times { class Distribution {
class BaseDistribution { public:
public: virtual ~Distribution() = default;
virtual ~BaseDistribution() = default;
BaseDistribution(std::string, double, std::mt19937_64&); Distribution(std::string, double, std::mt19937_64&);
explicit BaseDistribution(std::string, std::mt19937_64&); explicit Distribution(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; }
[[nodiscard]] std::string getName() const { return m_name; } [[nodiscard]] std::string getName() const { return m_name; };
virtual void setParameters(const std::unordered_map<std::string, double>&); 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() = 0;
virtual void draw_tau() = 0; virtual void draw_tau() = 0;
[[nodiscard]] double tau_wait() const; [[nodiscard]] double tau_wait() const;
[[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:
std::string m_name{"BaseDistribution"};
protected: double m_tau{1.};
std::string m_name{"BaseDistribution"}; double m_tau_jump{1.};
double m_tau{1.}; std::mt19937_64& m_rng;
double m_tau_jump{1.}; };
std::mt19937_64& m_rng;
};
}
#endif //RWSIM_TIMESBASE_H #endif //RWSIM_TIMESBASE_H

19
src/times/delta.cpp
View File

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

17
src/times/delta.h
View File

@ -3,16 +3,13 @@
#include "base.h" #include "base.h"
namespace times { class DeltaDistribution final : public Distribution {
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);
void initialize() override; void initialize() override;
void draw_tau() override; void draw_tau() override;
};
[[nodiscard]] std::string toString() const override;
};
}
#endif //RWSIM_TIMESDELTA_H #endif //RWSIM_TIMESDELTA_H

40
src/times/lognormal.cpp
View File

@ -1,31 +1,19 @@
#include "lognormal.h" #include "lognormal.h"
#include <cmath> #include <cmath>
namespace times { 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) {}
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) : Distribution(std::string("LogNormal"), rng) {}
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");
BaseDistribution::setParameters(parameters); Distribution::setParameters(parameters);
} }
std::unordered_map<std::string, double> LogNormalDistribution::getParameters() const { void LogNormalDistribution::initialize() {
auto parameter = BaseDistribution::getParameters(); m_distribution = std::lognormal_distribution(std::log(m_tau), m_sigma);
parameter["sigma"] = m_sigma; m_tau_jump = m_distribution(m_rng);
return parameter; }
}
void LogNormalDistribution::draw_tau() {
void LogNormalDistribution::initialize() { m_tau_jump = m_distribution(m_rng);
m_distribution = std::lognormal_distribution(std::log(m_tau), m_sigma);
m_tau_jump = m_distribution(m_rng);
}
void LogNormalDistribution::draw_tau() {
m_tau_jump = m_distribution(m_rng);
}
std::string LogNormalDistribution::toString() const {
return {"LogNormal/tau=" + std::to_string(m_tau) + "/sigma=" + std::to_string(m_sigma)};
}
} }

28
src/times/lognormal.h
View File

@ -3,25 +3,21 @@
#include "base.h" #include "base.h"
#include <random> #include <random>
#include <set>
namespace times { class LogNormalDistribution final : public Distribution {
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);
void setParameters(const std::unordered_map<std::string, double> &) 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; void initialize() override;
void draw_tau() override;
void initialize() override; private:
void draw_tau() override; double m_sigma{1};
std::lognormal_distribution<> m_distribution;
private: };
double m_sigma{1};
std::lognormal_distribution<> m_distribution;
};
}
#endif //LOGNORMAL_H #endif //LOGNORMAL_H

3
src/utils/functions.cpp
View File

@ -1,6 +1,9 @@
#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) {

52
src/utils/io.cpp
View File

@ -7,15 +7,11 @@
#include <complex> #include <complex>
#include <vector> #include <vector>
#include <iomanip> #include <iomanip>
#include <unordered_map>
#include <map> #include <map>
#include <string> #include <string>
#include <filesystem> #include <filesystem>
#include "../motions/base.h"
#include "../times/base.h"
class Motion;
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::string stripped_arg; std::string stripped_arg;
@ -25,7 +21,7 @@ std::pair<std::string, double> get_optional_parameter(std::vector<std::string>::
stripped_arg = it->substr(1, it->size()); stripped_arg = it->substr(1, it->size());
} }
std::transform(stripped_arg.begin(), stripped_arg.end(), stripped_arg.begin(), [](unsigned char c) { return std::tolower(c); }); std::transform(stripped_arg.begin(), stripped_arg.end(), stripped_arg.begin(), [](unsigned char c) { return std::tolower(c); });
const auto stripped_value = std::stod(*++it, nullptr); const auto stripped_value = std::stod(*(++it), nullptr);
return std::make_pair(stripped_arg, stripped_value); return std::make_pair(stripped_arg, stripped_value);
} }
@ -95,7 +91,7 @@ Arguments parse_args(const int argc, char* argv[]) {
std::unordered_map<std::string, double> read_parameter(const std::filesystem::path& infile) { std::unordered_map<std::string, double> read_parameter(const std::filesystem::path& infile) {
if (!exists(infile)) { if (!std::filesystem::exists(infile)) {
std::cerr << "File " << infile << " does not exist" << std::endl; std::cerr << "File " << infile << " does not exist" << std::endl;
exit(1); exit(1);
} }
@ -128,30 +124,16 @@ std::unordered_map<std::string, double> read_parameter(const std::filesystem::pa
} }
std::string make_directory(
const motions::BaseMotion& motion,
const times::BaseDistribution& distribution
) {
std::ostringstream path_name;
path_name << motion.toString() << "/" << distribution.toString();
if (!std::filesystem::create_directories(path_name.str())) {
std::cout << "Created directory " << path_name.str() << std::endl;
}
return path_name.str();
}
void save_parameter_to_file( void save_parameter_to_file(
const std::string& resulttype, const std::string& resulttype,
const std::string& directory, const std::string& motiontype,
const std::unordered_map<std::string, double>& parameter, const std::string& disttype,
const std::unordered_map<std::string, double>& optional std::unordered_map<std::string, double>& parameter,
std::unordered_map<std::string, double>& optional
) { ) {
std::ostringstream parameter_filename; std::ostringstream parameter_filename;
parameter_filename << directory << "/" << resulttype; parameter_filename << resulttype << "_" << motiontype << "_" << disttype;
parameter_filename << std::setprecision(6) << std::scientific; parameter_filename << std::setprecision(6) << std::scientific;
for (const auto& [key, value]: optional) { for (const auto& [key, value]: optional) {
@ -173,15 +155,15 @@ void save_parameter_to_file(
void save_data_to_file( void save_data_to_file(
const std::string& resulttype, const std::string& resulttype,
const std::string& directory, const std::string& motiontype,
const std::string& disttype,
const std::vector<double>& x, const std::vector<double>& x,
const std::map<double, std::vector<double>>& y, const std::map<double, std::vector<double>>& y,
const std::unordered_map<std::string, double>& optional std::unordered_map<std::string, double>& optional
) { ) {
// make file name // make file name
std::ostringstream datafile_name; std::ostringstream datafile_name;
datafile_name << directory << "/" << resulttype; datafile_name << resulttype << "_" << motiontype << "_" << disttype;
datafile_name << std::setprecision(6) << std::scientific; datafile_name << std::setprecision(6) << std::scientific;
for (const auto& [key, value]: optional) { for (const auto& [key, value]: optional) {
datafile_name << "_" << key << "=" << value; datafile_name << "_" << key << "=" << value;
@ -214,15 +196,15 @@ void save_data_to_file(
void save_data_to_file( void save_data_to_file(
const std::string& resulttype, const std::string& resulttype,
const std::string& directory, const std::string& motiontype,
const std::string& disttype,
const std::vector<double>& x, const std::vector<double>& x,
const std::vector<double>& y, const std::vector<double>& y,
const std::unordered_map<std::string, double>& optional std::unordered_map<std::string, double>& optional
) { ) {
// make file name // make file name
std::ostringstream datafile_name; std::ostringstream datafile_name;
datafile_name << directory << "/" << resulttype; datafile_name << resulttype << "_" << motiontype << "_" << disttype;
datafile_name << std::setprecision(6) << std::scientific; datafile_name << std::setprecision(6) << std::scientific;
for (const auto& [key, value]: optional) { for (const auto& [key, value]: optional) {
datafile_name << "_" << key << "=" << value; datafile_name << "_" << key << "=" << value;

11
src/utils/io.h
View File

@ -7,9 +7,6 @@
#include <filesystem> #include <filesystem>
#include <vector> #include <vector>
#include "../motions/base.h"
#include "../times/base.h"
struct Arguments { struct Arguments {
std::string parameter_file{}; std::string parameter_file{};
bool ste = false; bool ste = false;
@ -21,12 +18,12 @@ struct Arguments {
Arguments parse_args(int argc, char* argv[]); 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 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::string&, std::unordered_map<std::string, double>&, 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::string&, const std::vector<double>&, const std::map<double, std::vector<double>>&, std::unordered_map<std::string, double>&);
void save_data_to_file(const std::string&, const std::string&, const std::vector<double>&, const std::vector<double>&, const std::unordered_map<std::string, double>&); void save_data_to_file(const std::string&, const std::string&, const std::string&, const std::vector<double>&, const std::vector<double>&, std::unordered_map<std::string, double>&);
#endif #endif

16
start_sims.py
View File

@ -1,16 +0,0 @@
import numpy as np
from python.helpers import *
#Simulation parameter
motion = 'IsotropicAngle'
distribution = 'Delta'
parameter = {
"angle": [10, 109.47],
}
parameter = prepare_rw_parameter(parameter)
for variation in parameter:
print(f"\nRun RW for {motion}/{distribution} with arguments {variation}\n")
run_sims(motion, distribution, ste=True, spectrum=False, **variation)