diff --git a/.gitignore b/.gitignore index 4e122a1..ed13ab5 100644 --- a/.gitignore +++ b/.gitignore @@ -280,4 +280,14 @@ src/project2/buckling_beam_solver src/project2/off_diagonal_finder src/project2/scaling_tester -src/project2/*.txt \ No newline at end of file +src/project2/*.txt + +# Project 3 executables +src/project3/two_particles +src/project3/many_particles +src/project3/results/*.dat +src/project3/results/*/*.dat +src/project3/results/*/*/*.red +*.tar.gz +*_timing.txt +run_stats.json \ No newline at end of file diff --git a/projects/project2/main.pdf b/projects/project2/main.pdf index 41f6c95..23a1dc4 100644 Binary files a/projects/project2/main.pdf and b/projects/project2/main.pdf differ diff --git a/projects/project3/chapters/conclusion.tex b/projects/project3/chapters/conclusion.tex new file mode 100644 index 0000000..e69de29 diff --git a/projects/project3/chapters/introduction.tex b/projects/project3/chapters/introduction.tex new file mode 100644 index 0000000..2efb54a --- /dev/null +++ b/projects/project3/chapters/introduction.tex @@ -0,0 +1,5 @@ +Penning traps are useful appliances to capture charged particles without contact to any matter surrounding the particles \cite[11]{vogel_particle_2024}. This property of Penning traps is especially useful in the study of antimatter, as any contact to matter would lead to the direct annihilation of the antimatter. Therefore, Penning traps have become an essential tool in the antimatter research. Experiments like ALPHA-g \cite{bertsche_prospects_2018} and AEgIS \cite{scampoli_aegis_2014} at CERN have successfully deployed Penning traps for such applications. + +To plan future antimatter experiments it is essential to be able to simulate Penning traps. Without the ability of the behavior of particles within such a trap, a precise planning of the Penning trap and antimatter experiment is not possible. Key criterion for a successful simulation in this context is accuracy of the numerical simulation in terms of precise trajectories and key statistical properties of the captured particles, e.g. the kinetic energy. A second metric that is very important to the usefulness of such numerical simulations is the computational performance of the simulation as this decides whether the simulation is even possible or not with finite computing resources. To find the optimal balance between the finite computing power available and an accurate result, different integration algorithms will be compared to solve the ordinary differential equations at play in the Penning Trap. + +In \cref{sec:methods} the methods employed in the simulation will be expanded on. In detail the implementation of the source code will be discussed as well as the specific Penning trap used for the study, introduced. Thirdly the algorithms behind the numerical integration methods will be explained. \Cref{sec:results_and_discussion} will show the results obtained from our implementation and discuss their practical implications. This section will be split into four parts expanding on the accuracy of the trajectories, the energy conservation in the system, the computational performance of the different algorithms as well as a discussion on the limit of many particles in the trap. Lastly in \cref{sec:conclusion} our findings will be summarized and consequences of the results derived. \ No newline at end of file diff --git a/projects/project3/chapters/methods.tex b/projects/project3/chapters/methods.tex new file mode 100644 index 0000000..e47f9eb --- /dev/null +++ b/projects/project3/chapters/methods.tex @@ -0,0 +1,101 @@ +\subsection*{Penning trap} +Penning traps use a combination of a magnetic and electric field to confine the trajectory of a charged particle to a finite volume over long periods of time\cite[p.~9-10]{vogel_particle_2024}. Finite volume in this context refers to volumes in the order of \unit{\mm \cubed} up to \unit{\cm \cubed}. Long time ranges imply a period that exceeds the period of motion of the particle by multiple orders. For the studied case we will use an ideal Penning trap, i.e. there are no inhomogenities in the fields. The magnetic field will be defined as +\begin{equation} \label{eq:B_field} + B(\vec r, t) = B_0 \Theta\left(d - |\vec r |\right) \cdot \vec e_z, +\end{equation} +with the Heaviside function $\Theta(x)$, and we define two electric fields via the electric potentials +\begin{equation} \label{eq:V_stat} + \vstat(\vec r, t) = \frac{V_0 \cdot (2z^2 - x^2 - y^2) \Theta\left(d - |\vec r |\right)}{2d^2} +\end{equation} +and +\begin{equation} + \vdyn(\vec r, t) = (1 + f \cos \omega_V t) \cdot \vstat(\vec r, t). +\end{equation} + +\subsection*{Equations of motion} +Using Newton's equations of motion, $\ddot{\vec r} = \vec F / m$, we can derive ordinary differential equations for the position of a particle which is subject to \cref{eq:B_field,eq:V_stat}. As the magnetic field is parallel to the $z$-axis, the differential equation for the vertical component is straightforward to derive. Using the electric field resulting from $\vstat$ (see \cref{app:efield_equations} for details on all electric field components) we find +\begin{equation} + \ddot z + \omega_z^2 z \equiv \ddot z + \frac{2 q V_0}{m d^2} z = 0. +\end{equation} +In the transversal plane, the equations of motion are coupled due to the Lorentz force $\vec F = q \vec v \times \vec B$. Using the electric field resulting from $\vstat$ (see \cref{app:efield_equations} for details on all electric field components) we find +\begin{align} + \ddot x - \omega_0 \dot y - \frac{1}{2} \omega_z^2 x &\equiv \ddot x - \frac{q B_0}{m} \dot y - \frac{q V_0}{m d^2} x = 0, \label{eq:transv_eom_x}\\ + \ddot y + \omega_0 \dot x - \frac{1}{2} \omega_z^2 y &\equiv \ddot y + \frac{q B_0}{m} \dot x - \frac{q V_0}{m d^2} y = 0 \label{eq:transv_eom_y}. +\end{align} +Using the sum of equations: $\eqref{eq:transv_eom_x} + i \eqref{eq:transv_eom_y}$, we can use the general definition $f(t) \equiv x(t) + i y(t)$ to derive a single complex differential equation for the transversal motion +\begin{equation} \label{eq:complex_transv_eom} + \ddot f + i \omega_0 \dot f - \frac{1}{2} \omega_z^2 f = 0. +\end{equation} +After solving the differential equation, the real and imaginary parts of $f(t)$ will correspond to the $x$ and $y$ components of the trajectory, respectively. +\Cref{eq:complex_transv_eom} is a damped harmonic oscillator equation with the general solution +\begin{equation} \label{eq:general_transv_solution} + f(t) = A_+ e^{-i (\omega_+ t + \phi_+)} + A_- e^{-i (\omega_- t + \phi_-)}, +\end{equation} +with the two characteristic frequencies +\begin{equation} + \omega_\pm = \frac{\omega_0 \pm \sqrt{\omega_0^2 - 2 \omega_z^2}}{2}. +\end{equation} +This solution is valid for $\omega_0^2 > 2 \omega_z^2$. In the case of $\omega_0^2 \leq \omega_z^2$, the oscillation frequencies become complex, leading to exponential growth of the trajectory in the transversal plane. This condition is known as the stability criterion for a Penning trap \cite[p.~67]{vogel_particle_2024}. In terms of the physical parameters of the trap, the stability criterion can be expressed as +\begin{equation} + \frac{q}{m} B_0^2 > \frac{2 V_0}{d^2}. +\end{equation} +Given this condition, the solution for the transversal motion will be bound, i.e. $|f(t)| < \infty$ for all $t > 0$. Given the stability criterion, the two terms in \cref{eq:general_transv_solution} can be identified as sinusoidal motions with amplitudes $A_\pm$. If the two terms are in phase the upper limit of $|f| \equiv R_+$ will be $R_+ = A_+ + A_-$, while the lower limit for the case of antiparallel phases will be the absolute of the difference in the two amplitudes, $R_- = |A_+ - A_-|$. A analytical solution for a special case of initial conditions is derived in \cref{app:special_case_analytical_solution}. + +\subsubsection*{Multiple particles} +For the case of multiple particles, it is important to include the particle-particle interactions. For this analysis we will only consider the Coulomb interaction between the particles, neglecting magnetic interactions. The force on particle $i$ due to all other particles $j$ is then given by +\begin{equation} + \vec F_i = \frac{q_i}{4 \pi \epsilon_0} \sum_{j \neq i} \frac{q_j (\vec r_j - \vec r_i)}{|\vec r_j - \vec r_i|^3}. +\end{equation} + +\subsection*{Numerical Integration Methods} +\subsubsection*{Forward Euler method} +The simplest numerical method to solve the ordinary differential equation system is the Forward Euler method. It is a first order method, meaning that the local truncation error per step is on the order of $\mathcal{O}(h^2)$, with $h$ being the step size. The global error after $N$ steps is therefore on the order of $\mathcal{O}(h)$. The method is explicit, meaning that the state of the system at the next time step can be calculated directly from the current state. Given a general ordinary differential equation of the form +\begin{equation} + \dot y(t) = f(t, y(t)), +\end{equation} +the Forward Euler method updates the state of the system as follows: +\begin{equation} + y_{n+1} = y_n + h f(t_n, y_n). +\end{equation} +For our second order differential equations, we first rewrite them as a system of first order equations. The next state of the system is then calculated using the current position and velocity as +\begin{align} + \vec r_{n+1} &= \vec r_n + h \vec v_n, \\ + \vec v_{n+1} &= \vec v_n + h \frac{\vec F(\vec r_n, \vec v_n, t_n)}{m}. +\end{align} +This method is computationally inexpensive, as we only need a single force evaluation per time step. However, it is known to be not very accurate. +\subsubsection*{Runge-Kutta 4 method} +The Runge-Kutta 4 (RK4) method is a popular and more accurate method for solving ordinary differential equations. It is a fourth-order method, meaning that the local truncation error per step is on the order of $\mathcal{O}(h^5)$, and the global error after $N$ steps is on the order of $\mathcal{O}(h^4)$. The RK4 method calculates the next state of the system using a weighted average of four different estimates of the slope (the derivative) at different points within the time step. Given a general ordinary differential equation of the form +\begin{equation} + \dot y(t) = f(t, y(t)), +\end{equation} +the RK4 method updates the state of the system as follows: +\begin{align} + k_1 &= h f(t_n, y_n), \\ + k_2 &= h f\left(t_n + \frac{h}{2}, y_n + \frac{k_1}{2}\right), \\ + k_3 &= h f\left(t_n + \frac{h}{2}, y_n + \frac{k_2}{2}\right), \\ + k_4 &= h f(t_n + h, y_n + k_3), \\ + y_{n+1} &= y_n + \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4). +\end{align} +As this method is again a generalized solver for ODE of the first order, we again rewrite our second order differential equations as a system of first order equations. + +\subsubsection*{Velocity-Verlet method} +In contrast to the previous two methods, the Velocity-Verlet method is a symplectic integrator which is specifically designed for second order differential equations of the form +\begin{equation} + \ddot{\vec r}(t) = \frac{\vec F(\vec r(t), t)}{m}. +\end{equation} +The method is time-reversible and conserves energy better over long time periods compared to non-symplectic methods like Forward Euler and RK4. The Velocity-Verlet method updates the position and velocity of the system as follows: +\begin{align} + \vec r_{n+1} &= \vec r_n + h \vec v_n + \frac{h^2}{2} \vec a_n, \\ + \vec a_{n+1} &= \frac{\vec F(\vec r_{n+1}, t_{n+1})}{m}, \\ + \vec v_{n+1} &= \vec v_n + \frac{h}{2} (\vec a_n + \vec a_{n+1}), +\end{align} +It is important to note that the Velocity-Verlet integrator requires forces to be independent of velocity. In the case, that the algorithm is applicable, it is a very efficient method, as it only requires a single force evaluation per time step, while still being a second order method with a local truncation error per step on the order of $\mathcal{O}(h^3)$ and a global error after $N$ steps on the order of $\mathcal{O}(h^2)$. + +\subsection*{Code Structure} +The basic framework for the numerical analysis is based on a \texttt{PenningTrap} class, which contains all the particles present in the trap, as well as a parametrization of the electric and magnetic fields. The particles are represented by a \texttt{Particle} class, which contains the physical properties of the particle, as well as its current position and velocity. With this information, the \texttt{PenningTrap} class can calculate the forces acting on each particle, including the external fields and the particle-particle interactions. The particle-particle interactions can be toggled on and off, allowing for a simulation of both scenarios. The external fields can be modified by supplying a field-method of the form \texttt{external\_field(const arma::vec\& r, double t, const PenningTrap\& trap)}. The reference to the \texttt{PenningTrap} allows for the parameters of the field to be stored in the trap object. The implementation of the \texttt{Particle} and \texttt{PenningTrap} classes can be found in \texttt{/src/project3/include/classes.hpp} of the project repository, as well as in the corresponding source file \texttt{/src/project3/src/classes.cpp}. + +\subsection*{Numerical Methods Implementation} +All numerical methods are implemented as classes inheriting from a base class \texttt{Solver}. The base class contains a reference to the \texttt{PenningTrap} object, as well as the time step size. The recording of particle properties over time, like postion and velocity is part of the general \texttt{Solver} class. Each derived class implements the \texttt{step()} method, which updates the state of the system by one time step using the respective numerical method. The implementation of the \texttt{Solver} class and its derived classes can be found in \texttt{/src/project3/include/solvers.hpp} of the project repository, as well as in the corresponding source file \texttt{/src/project3/src/solvers.cpp}. The following solvers are implemented:\texttt{Forward\-Euler\-Solver}, \texttt{RK4\-Solver}, \texttt{Velocity\-Verlet\-Solver} and \texttt{Analytical\-Solver}, which implements the analytical solution for a special case of initial conditions (see \cref{app:special_case_analytical_solution}). + +\subsection*{Tools} +\textcolor{red}{TODO} \ No newline at end of file diff --git a/projects/project3/chapters/results.tex b/projects/project3/chapters/results.tex new file mode 100644 index 0000000..7ea86a7 --- /dev/null +++ b/projects/project3/chapters/results.tex @@ -0,0 +1,9 @@ +\subsection{Numerical Accuracy of Trajectories} + +\subsection{Energy Conservation} + + +\subsection{Performance and Efficiency} + +\subsection{Many-Body Simulations} + diff --git a/projects/project3/main.pdf b/projects/project3/main.pdf new file mode 100644 index 0000000..f84eb78 Binary files /dev/null and b/projects/project3/main.pdf differ diff --git a/projects/project3/main.tex b/projects/project3/main.tex new file mode 100644 index 0000000..128d958 --- /dev/null +++ b/projects/project3/main.tex @@ -0,0 +1,156 @@ +% USEFUL LINKS: +% ------------- +% +% - UiO LaTeX guides: https://www.mn.uio.no/ifi/tjenester/it/hjelp/latex/ +% - Mathematics: https://en.wikibooks.org/wiki/LaTeX/Mathematics +% - Physics: https://ctan.uib.no/macros/latex/contrib/physics/physics.pdf +% - Basics of Tikz: https://en.wikibooks.org/wiki/LaTeX/PGF/Tikz +% - All the colors! https://en.wikibooks.org/wiki/LaTeX/Colors +% - How to make tables: https://en.wikibooks.org/wiki/LaTeX/Tables +% - Code listing styles: https://en.wikibooks.org/wiki/LaTeX/Source_Code_Listings +% - \includegraphics https://en.wikibooks.org/wiki/LaTeX/Importing_Graphics +% - Learn more about figures: https://en.wikibooks.org/wiki/LaTeX/Floats,_Figures_and_Captions +% - Automagic bibliography: https://en.wikibooks.org/wiki/LaTeX/Bibliography_Management (this one is kinda difficult the first time) +% +% (This document is of class "revtex4-1", the REVTeX Guide explains how the class works) +% REVTeX Guide: http://www.physics.csbsju.edu/370/papers/Journal_Style_Manuals/auguide4-1.pdf +% +% +% COMPILING THE .pdf FILE IN THE LINUX TERMINAL +% --------------------------------------------- +% +% [terminal]$ pdflatex report_example.tex +% +% Run the command twice, always. +% +% When using references, footnotes, etc. you should run the following chain of commands: +% +% [terminal]$ pdflatex report_example.tex +% [terminal]$ bibtex report_example +% [terminal]$ pdflatex report_example.tex +% [terminal]$ pdflatex report_example.tex +% +% This series of commands can of course be gathered into a single-line command: +% [terminal]$ pdflatex report_example.tex && bibtex report_example.aux && pdflatex report_example.tex && pdflatex report_example.tex +% +% ---------------------------------------------------- + + + +% \documentclass[english,notitlepage,reprint,nofootinbib]{revtex4-2} % defines the basic parameters of the document +\documentclass[english,notitlepage,reprint,nofootinbib]{revtex4-2} % defines the basic parameters of the document +% If you want a single-column, remove "reprint" +\usepackage{silence} +\WarningFilter{revtex4-2}{Repair the float} + +% Allows special characters (including æøå) +\usepackage[utf8]{inputenc} +\usepackage[english]{babel} + +% Note that you may need to download some of these packages manually, it depends on your setup. +% It may be usefult to download TeXMaker, because it includes a large library of the most common packages. + +\usepackage{amssymb} % mathematical symbols (physics imports amsmath) +\usepackage{amsmath} +\usepackage{graphicx} % include graphics such as plots +\usepackage{xcolor} % set colors +\usepackage{hyperref} % automagic cross-referencing +\usepackage{listings} % display code +\usepackage{subfigure} % imports a lot of cool and useful figure commands +% \usepackage{float} +%\usepackage[section]{placeins} +\usepackage{algorithm} +\usepackage[noend]{algpseudocode} +\usepackage{subfigure} +\usepackage{tikz} +\usepackage{cleveref} +\usepackage{siunitx} +\usetikzlibrary{quantikz} +% defines the color of hyperref objects +% Blending two colors: blue!80!black = 80% blue and 20% black +\hypersetup{ % this is just my personal choice, feel free to change things + colorlinks, + linkcolor={red!50!black}, + citecolor={blue!50!black}, + urlcolor={blue!80!black}} + + +% =========================================== +\newcommand{\vstat}{V_\mathrm{stat.}} +\newcommand{\vdyn}{V_\mathrm{dyn.}} +\newcommand{\estat}{E_\mathrm{stat.}} +\newcommand{\edyn}{E_\mathrm{dyn.}} +\newcommand{\estatc}[1]{E_{\mathrm{stat.},#1}} +\newcommand{\edync}[1]{E_{\mathrm{dyn.},#1}} + +\begin{document} + +\title{Accuracy and Stability of Numerical Integration Methods in Penning Trap Simulations} % self-explanatory +\author{Lars Bogner} % self-explanatory +\date{\today} % self-explanatory +\noaffiliation % ignore this, but keep it. + +%This is how we create an abstract section. +\begin{abstract} + \textcolor{red}{TODO: ABSTRACT HERE.} +\end{abstract} +\maketitle + + +% =========================================== +\section{Introduction} +\input{chapters/introduction} + +% =========================================== +\section{Methods}\label{sec:methods} +% +\input{chapters/methods} + + +% =========================================== +\section{Results and discussion}\label{sec:results_and_discussion} +% +\input{chapters/results} + +% =========================================== +\section{Conclusion}\label{sec:conclusion} +% +\input{chapters/conclusion} + +% =========================================== +\appendix +\section{Electric Field Equations} \label{app:efield_equations} +Here we derive the electric field equations used in \cref{sec:methods}. The electric field following from $\vstat$ defined in \cref{eq:V_stat} is given by +\begin{align} + \estatc{x} &= -\frac{\partial \vstat}{\partial x} = \frac{V_0}{d^2} x, \\ + \estatc{y} &= -\frac{\partial \vstat}{\partial y} = \frac{V_0}{d^2} y, \\ + \estatc{z} &= -\frac{\partial \vstat}{\partial z} = -\frac{2 V_0}{d^2} z, +\end{align} +for all $|\vec r| < d$ and zero otherwise. The electric field following from $\vdyn$ is simply +\begin{equation} + \edyn = (1 + f \cos \omega_V t) \cdot \estat. +\end{equation} + +\section{Special Case Analytical Solution} \label{app:special_case_analytical_solution} +In the special case of a single particle in the potential $\vstat$, with the inital conditions $\vec{r}(t_0) = (x_0, 0, z_0)$ and $\vec{v}(t_0) = (0, v_0, 0)$, we can derive an analytical solution for the trajectory. The solution for the $z$ component is straightforward +\begin{equation} + z(t) = z_0 \cos(\omega_z t). +\end{equation} +From the equation system +\begin{align} + f(0) = A_+ e^{-i \phi_+} + A_- e^{-i \phi_-} &\equiv x_0, \\ + \dot f(0) = -i \omega_+ A_+ e^{-i \phi_+} - i \omega_- A_- e^{-i \phi_-} &\equiv i v_0, +\end{align} +we derive $\phi_\pm = 0$ as $f(0) \in \mathbb{R}$ and therefore $A_+ + A_- = x_0$ and $\omega_+ A_+ + \omega_- A_- = -v_0$. Solving this system of equations for $A_\pm$ we find +\begin{equation} + A_\pm =\pm \frac{v_0 + \omega_\mp x_0}{\omega_- - \omega_+}. +\end{equation} + +% =========================================== +\onecolumngrid +% \bibliographystyle{apalike} +\bibliographystyle{unsrt} +\bibliography{ref} + + +\end{document} \ No newline at end of file diff --git a/projects/project3/mainNotes.bib b/projects/project3/mainNotes.bib new file mode 100644 index 0000000..e69de29 diff --git a/projects/project3/ref.bib b/projects/project3/ref.bib new file mode 100644 index 0000000..299e0f8 --- /dev/null +++ b/projects/project3/ref.bib @@ -0,0 +1,187 @@ + +@article{harris_array_2020, + title = {Array programming with {NumPy}}, + volume = {585}, + url = {https://doi.org/10.1038/s41586-020-2649-2}, + doi = {10.1038/s41586-020-2649-2}, + pages = {357--362}, + number = {7825}, + journaltitle = {Nature}, + author = {Harris, Charles R. and Millman, K. Jarrod and Walt, Stéfan J. van der and Gommers, Ralf and Virtanen, Pauli and Cournapeau, David and Wieser, Eric and Taylor, Julian and Berg, Sebastian and Smith, Nathaniel J. and Kern, Robert and Picus, Matti and Hoyer, Stephan and Kerkwijk, Marten H. van and Brett, Matthew and Haldane, Allan and Río, Jaime Fernández del and Wiebe, Mark and Peterson, Pearu and Gérard-Marchant, Pierre and Sheppard, Kevin and Reddy, Tyler and Weckesser, Warren and Abbasi, Hameer and Gohlke, Christoph and Oliphant, Travis E.}, + date = {2020-09}, + note = {Publisher: Springer Science and Business Media {LLC}}, +} + +@article{hunter_matplotlib_2007, + title = {Matplotlib: A 2D graphics environment}, + volume = {9}, + doi = {10.1109/MCSE.2007.55}, + abstract = {Matplotlib is a 2D graphics package used for Python for application development, interactive scripting, and publication-quality image generation across user interfaces and operating systems.}, + pages = {90--95}, + number = {3}, + journaltitle = {Computing in Science \& Engineering}, + author = {Hunter, J. D.}, + date = {2007}, + note = {Publisher: {IEEE} {COMPUTER} {SOC}}, +} + +@article{pedregosa_scikit-learn_2011, + title = {Scikit-learn: Machine Learning in Python}, + volume = {12}, + pages = {2825--2830}, + journaltitle = {Journal of Machine Learning Research}, + author = {Pedregosa, F. and Varoquaux, G. and Gramfort, A. and Michel, V. and Thirion, B. and Grisel, O. and Blondel, M. and Prettenhofer, P. and Weiss, R. and Dubourg, V. and Vanderplas, J. and Passos, A. and Cournapeau, D. and Brucher, M. and Perrot, M. and Duchesnay, E.}, + date = {2011}, +} + +@inreference{noauthor_stochastic_2025, + title = {Stochastic gradient descent}, + rights = {Creative Commons Attribution-{ShareAlike} License}, + url = {https://en.wikipedia.org/w/index.php?title=Stochastic_gradient_descent&oldid=1309164477}, + abstract = {Stochastic gradient descent (often abbreviated {SGD}) is an iterative method for optimizing an objective function with suitable smoothness properties (e.g. differentiable or subdifferentiable). It can be regarded as a stochastic approximation of gradient descent optimization, since it replaces the actual gradient (calculated from the entire data set) by an estimate thereof (calculated from a randomly selected subset of the data). Especially in high-dimensional optimization problems this reduces the very high computational burden, achieving faster iterations in exchange for a lower convergence rate. +The basic idea behind stochastic approximation can be traced back to the Robbins–Monro algorithm of the 1950s. Today, stochastic gradient descent has become an important optimization method in machine learning.}, + booktitle = {Wikipedia}, + urldate = {2025-09-22}, + date = {2025-09-02}, + langid = {english}, + note = {Page Version {ID}: 1309164477}, + file = {Snapshot:/home/lars/Zotero/storage/B6CVR59B/index.html:text/html}, +} + +@unpublished{elstner_lecture_2025, + location = {Karlsruhe Institute for Technology, Karlsruhe}, + title = {Lecture: Machine Learning for Chemistry}, + type = {Lecture}, + howpublished = {Lecture}, + author = {Elstner, Marcus and Kubar, Tomas}, + date = {2025-05-20}, + langid = {german}, + file = {PDF:/home/lars/Zotero/storage/5LSLJMK8/Elstner and Kubar - 2025 - Lecture Machine Learning for Chemistry.pdf:application/pdf}, +} + +@online{lekhansh_lasso_2024, + title = {Lasso vs. Ridge Regression: A Detailed Comparison}, + url = {https://medium.com/@tyagi.lekhansh/lasso-vs-ridge-regression-a-detailed-comparison-140f7832c624}, + shorttitle = {Lasso vs. Ridge Regression}, + abstract = {In the realm of regression analysis, Lasso (Least Absolute Shrinkage and Selection Operator) and Ridge Regression are two popular…}, + titleaddon = {Medium}, + author = {Lekhansh}, + urldate = {2025-09-22}, + date = {2024-09-04}, + langid = {english}, +} + +@book{hastie_elements_2009, + location = {New York, {NY}}, + title = {The Elements of Statistical Learning}, + rights = {http://www.springer.com/tdm}, + isbn = {978-0-387-84857-0 978-0-387-84858-7}, + url = {http://link.springer.com/10.1007/978-0-387-84858-7}, + series = {Springer Series in Statistics}, + publisher = {Springer}, + author = {Hastie, Trevor and Tibshirani, Robert and Friedman, Jerome}, + urldate = {2025-09-22}, + date = {2009}, + doi = {10.1007/978-0-387-84858-7}, + keywords = {Averaging, Boosting, classification, clustering, data mining, machine learning, Projection pursuit, Random Forest, supervised learning, Support Vector Machine, unsupervised learning}, + file = {Full Text PDF:/home/lars/Zotero/storage/D3N4DVY9/Hastie et al. - 2009 - The Elements of Statistical Learning.pdf:application/pdf}, +} + +@book{goodfellow_deep_2016, + title = {Deep Learning}, + publisher = {{MIT} Press}, + author = {Goodfellow, Ian and Bengio, Yoshua and Courville, Aaron}, + date = {2016}, +} + +@book{bishop_pattern_2006, + location = {New York}, + title = {Pattern recognition and machine learning}, + isbn = {978-0-387-31073-2}, + series = {Information science and statistics}, + publisher = {Springer}, + author = {Bishop, Christopher M.}, + date = {2006}, + langid = {english}, + file = {PDF:/home/lars/Zotero/storage/9H5W9BGC/Bishop - 2006 - Pattern recognition and machine learning.pdf:application/pdf}, +} + +@online{roskam_supercritical_2025, + title = {Supercritical Fluid Chromatography for Chiral Analysis, Part 1: Theoretical Background {\textbar} {LCGC} International}, + url = {https://www.chromatographyonline.com/view/supercritical-fluid-chromatography-for-chiral-analysis-part-1-theoretical-background}, + shorttitle = {Supercritical Fluid Chromatography for Chiral Analysis, Part 1}, + abstract = {With the substantial developments carried out over the past years in instrumentation, columns, and detector hyphenation, the interest in chiral supercritical fluid chromatography ({SFC}) has been steadily growing in various fields. In the first part of this review article, the theoretical advantages, technological developments, and common practices in chiral {SFC} are discussed.}, + author = {Roskam, Gerry and Velde, Bas van de and Gargano, Andrea and Kohler, Isabelle}, + urldate = {2025-10-05}, + date = {2025-10-05}, + langid = {english}, + file = {Snapshot:/home/lars/Zotero/storage/LWSHJMVJ/supercritical-fluid-chromatography-for-chiral-analysis-part-1-theoretical-background.html:text/html}, +} + +@article{furet_first_2022, + title = {The First Class of Small Molecules Potently Disrupting the {YAP}‐{TEAD} Interaction by Direct Competition}, + volume = {17}, + doi = {10.1002/cmdc.202200303}, + abstract = {Inhibition of the {YAP}‐{TEAD} protein‐protein interaction is an attractive therapeutic concept under intense investigation with the objective to treat cancers associated with a dysregulation of the Hippo pathway. However, owing to the very extended surface of interaction of the two proteins, the identification of small drug‐like molecules able to efficiently prevent {YAP} from binding to {TEAD} by direct competition has been elusive so far. We disclose here the discovery of the first class of small molecules potently inhibiting the {YAP}‐{TEAD} interaction by binding at one of the main interaction sites of {YAP} at the surface of {TEAD}. These inhibitors, providing a path forward to pharmacological intervention in the Hippo pathway, evolved from a weakly active virtual screening hit advanced to high potency by structure‐based design.}, + journaltitle = {{ChemMedChem}}, + shortjournal = {{ChemMedChem}}, + author = {Furet, Pascal and Bordas, Vincent and Douget, Mickaël and Salem, Bahaa and Mesrouze, Yannick and Imbach‐Weese, Patricia and Sellner, Holger and Vögtle, Markus and Soldermann, Nicolas and Chapeau, Emilie and Wartmann, Markus and Scheufler, Clemens and Fernández, César and Kallen, Joerg and Guagnano, Vito and Chene, Patrick and Schmelzle, Tobias}, + date = {2022-09-02}, + file = {Full Text PDF:/home/lars/Zotero/storage/J399MX84/Furet et al. - 2022 - The First Class of Small Molecules Potently Disrupting the YAP‐TEAD Interaction by Direct Competitio.pdf:application/pdf}, +} + +@book{vogel_particle_2024, + location = {Cham}, + title = {Particle Confinement in Penning Traps: An Introduction}, + volume = {126}, + rights = {https://www.springernature.com/gp/researchers/text-and-data-mining}, + isbn = {978-3-031-55419-3 978-3-031-55420-9}, + url = {https://link.springer.com/10.1007/978-3-031-55420-9}, + series = {Springer Series on Atomic, Optical, and Plasma Physics}, + shorttitle = {Particle Confinement in Penning Traps}, + publisher = {Springer International Publishing}, + author = {Vogel, Manuel}, + urldate = {2025-10-09}, + date = {2024}, + langid = {english}, + doi = {10.1007/978-3-031-55420-9}, + keywords = {Confined ions and plasmas, Highly charged ions, Ion trapping, Laser Cooling, Magnetic moments, Particle confinement, Penning traps, Precision spectroscopy, Resistive Cooling, Stored ions, Trapped charged particles, Unified notation Penning traps}, + file = {Full Text PDF:/home/lars/Zotero/storage/SI8SF7RF/Vogel - 2024 - Particle Confinement in Penning Traps An Introduction.pdf:application/pdf}, +} + +@article{scampoli_aegis_2014, + title = {The {AEgIS} experiment at {CERN} for the measurement of antihydrogen gravity acceleration}, + volume = {29}, + issn = {0217-7323}, + url = {https://www.worldscientific.com/doi/abs/10.1142/S0217732314300171}, + doi = {10.1142/S0217732314300171}, + abstract = {The Antihydrogen Experiment: Gravity, Interferometry, Spectroscopy ({AEgIS}) experiment is conducted by an international collaboration based at {CERN} whose aim is to perform the first direct measurement of the gravitational acceleration of antihydrogen in the local field of the Earth, with Δg/g = 1\% precision as a first achievement. The idea is to produce cold (100 {mK}) antihydrogen through a pulsed charge exchange reaction by overlapping clouds of antiprotons, from the Antiproton Decelerator ({AD}) and positronium atoms inside a Penning trap. The antihydrogen has to be produced in an excited Rydberg state to be subsequently accelerated to form a beam. The deflection of the antihydrogen beam can then be measured by using a moiré deflectometer coupled to a position sensitive detector to register the impact point of the anti-atoms through the vertex reconstruction of their annihilation products. After being approved in late 2008, {AEgIS} started taking data in a commissioning phase in 2012. This paper presents an outline of the experiment with a brief overview of its physics motivation and of the state-of-the-art of the g measurement on antimatter. Particular attention is given to the current status of the emulsion-based position detector needed to measure the sag in {AEgIS}.}, + pages = {1430017}, + number = {17}, + journaltitle = {Modern Physics Letters A}, + shortjournal = {Mod. Phys. Lett. A}, + author = {Scampoli, Paola and Storey, James}, + urldate = {2025-10-09}, + date = {2014-06-07}, + note = {Publisher: World Scientific Publishing Co.}, + keywords = {Antihydrogen, gravity, high resolution tracking detector}, +} + +@article{bertsche_prospects_2018, + title = {Prospects for comparison of matter and antimatter gravitation with {ALPHA}-g}, + volume = {376}, + url = {https://royalsocietypublishing.org/doi/full/10.1098/rsta.2017.0265}, + doi = {10.1098/rsta.2017.0265}, + abstract = {The {ALPHA} experiment has recently entered an expansion phase of its experimental programme, driven in part by the expected benefits of conducting experiments in the framework of the new {AD} + {ELENA} antiproton facility at {CERN}. With antihydrogen trapping now a routine operation in the {ALPHA} experiment, the collaboration is leading progress towards precision atomic measurements on trapped antihydrogen atoms, with the first excitation of the 1S–2S transition and the first measurement of the antihydrogen hyperfine spectrum (Ahmadi et al. 2017 Nature 541, 506–510 (doi:10.1038/nature21040); Nature 548, 66–69 (doi:10.1038/nature23446)). We are building on these successes to extend our physics programme to include a measurement of antimatter gravitation. We plan to expand a proof-of-principle method (Amole et al. 2013 Nat. Commun. 4, 1785 (doi:10.1038/ncomms2787)), first demonstrated in the original {ALPHA} apparatus, and perform a precise measurement of antimatter gravitational acceleration with the aim of achieving a test of the weak equivalence principle at the 1\% level. The design of this apparatus has drawn from a growing body of experience on the simulation and verification of antihydrogen orbits confined within magnetic-minimum atom traps. The new experiment, {ALPHA}-g, will be an additional atom-trapping apparatus located at the {ALPHA} experiment with the intention of measuring antihydrogen gravitation. + +This article is part of the Theo Murphy meeting issue ‘Antiproton physics in the {ELENA} era’.}, + pages = {20170265}, + number = {2116}, + journaltitle = {Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences}, + author = {Bertsche, W. A.}, + urldate = {2025-10-09}, + date = {2018-02-19}, + note = {Publisher: Royal Society}, + keywords = {antigravity, antihydrogen, antimatter, {CPT}, gravity, Lorentz invariance}, + file = {Full Text PDF:/home/lars/Zotero/storage/7YM83XXD/Bertsche - 2018 - Prospects for comparison of matter and antimatter gravitation with ALPHA-g.pdf:application/pdf}, +} diff --git a/pyproject.toml b/pyproject.toml index c101e28..f4bfa6b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -8,4 +8,5 @@ dependencies = [ "numpy>=2.3.2", "pandas>=2.3.2", "scipy>=1.16.1", + "typer>=0.19.2", ] diff --git a/src/project3/generate_results.sh b/src/project3/generate_results.sh new file mode 100755 index 0000000..4821636 --- /dev/null +++ b/src/project3/generate_results.sh @@ -0,0 +1,87 @@ +#!/bin/bash + +PUEUE=true +NORMAL_MANY_PARTICLE_SIMULATION=true +OSCILLATING_FIELD_SIMULATION=false + +echo "Removing old results..." +mkdir -p results +rm results/*.dat + +echo "Compiling binaries..." +rm -f two_particles many_particles +g++ -O3 -larmadillo -I ./include src/* two_particles.cpp -o two_particles +g++ -O3 -larmadillo -I ./include src/* many_particles.cpp -o many_particles +echo "Finished compiling." + +if $PUEUE; then + echo "Using pueue to generate results in parallel..." + systemctl --user start pueued + pueue parallel 8 + pueue clean + COMMAND_PREFIX="pueue add " +else + echo "Generating results sequentially..." + COMMAND_PREFIX="" +fi + +for N in 4000 8000 16000 32000 +do + for interactions in "" "--disable-interactions" + do + for n in 1 2 + do + $COMMAND_PREFIX ./two_particles -N $N -n $n $interactions + $COMMAND_PREFIX ./two_particles -N $N -n $n $interactions -E + $COMMAND_PREFIX ./two_particles -N $N -n $n $interactions -V + done + done + # The analytical solution is only valid for one particle (second particle has v_z) without interactions + $COMMAND_PREFIX ./two_particles -N $N -n 1 --disable-interactions -A +done + +if $NORMAL_MANY_PARTICLE_SIMULATION; then + for N in 4000 8000 16000 32000 + do + for interactions in "" "--disable-interactions" + do + for n in 10 50 100 + do + $COMMAND_PREFIX ./many_particles -N $N -n $n $interactions + $COMMAND_PREFIX ./many_particles -N $N -n $n $interactions -E + $COMMAND_PREFIX ./many_particles -N $N -n $n $interactions -V + done + done + done +fi + +if $OSCILLATING_FIELD_SIMULATION; then + N=40000 + n=100 + reduced_output="--reduced-output" + interactions="--disable-interactions" + for amplitude in 0.1 0.4 0.7 + do + for frequency in $(seq 0.2 0.02 2.5) + do + $COMMAND_PREFIX ./many_particles -N $N -t -n $n $interactions -O -f $amplitude -w $frequency $reduced_output + $COMMAND_PREFIX ./many_particles -N $N -t -n $n $interactions -V -O -f $amplitude -w $frequency $reduced_output + done + done + interactions="" + for amplitude in 0.1 0.4 0.7 + do + for frequency in $(seq 1.0 0.005 1.8) + do + $COMMAND_PREFIX ./many_particles -N $N -t -n $n $interactions -O -f $amplitude -w $frequency $reduced_output + $COMMAND_PREFIX ./many_particles -N $N -t -n $n $interactions -V -O -f $amplitude -w $frequency $reduced_output + done + done +fi + + +echo "All tasks added to pueue." + +echo "Use 'pueue status' to check the status of the tasks." + +pueue status \ No newline at end of file diff --git a/src/project3/include/argparse/argparse.hpp b/src/project3/include/argparse/argparse.hpp new file mode 100644 index 0000000..06d30fd --- /dev/null +++ b/src/project3/include/argparse/argparse.hpp @@ -0,0 +1,2589 @@ +/* + __ _ _ __ __ _ _ __ __ _ _ __ ___ ___ + / _` | '__/ _` | '_ \ / _` | '__/ __|/ _ \ Argument Parser for Modern C++ +| (_| | | | (_| | |_) | (_| | | \__ \ __/ http://github.com/p-ranav/argparse + \__,_|_| \__, | .__/ \__,_|_| |___/\___| + |___/|_| + +Licensed under the MIT License . +SPDX-License-Identifier: MIT +Copyright (c) 2019-2022 Pranav Srinivas Kumar +and other contributors. + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. +*/ +#pragma once + +#include + +#ifndef ARGPARSE_MODULE_USE_STD_MODULE +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#endif + +#ifndef ARGPARSE_CUSTOM_STRTOF +#define ARGPARSE_CUSTOM_STRTOF strtof +#endif + +#ifndef ARGPARSE_CUSTOM_STRTOD +#define ARGPARSE_CUSTOM_STRTOD strtod +#endif + +#ifndef ARGPARSE_CUSTOM_STRTOLD +#define ARGPARSE_CUSTOM_STRTOLD strtold +#endif + +namespace argparse { + +namespace details { // namespace for helper methods + +template +struct HasContainerTraits : std::false_type {}; + +template <> struct HasContainerTraits : std::false_type {}; + +template <> struct HasContainerTraits : std::false_type {}; + +template +struct HasContainerTraits< + T, std::void_t().begin()), + decltype(std::declval().end()), + decltype(std::declval().size())>> : std::true_type {}; + +template +inline constexpr bool IsContainer = HasContainerTraits::value; + +template +struct HasStreamableTraits : std::false_type {}; + +template +struct HasStreamableTraits< + T, + std::void_t() << std::declval())>> + : std::true_type {}; + +template +inline constexpr bool IsStreamable = HasStreamableTraits::value; + +constexpr std::size_t repr_max_container_size = 5; + +template std::string repr(T const &val) { + if constexpr (std::is_same_v) { + return val ? "true" : "false"; + } else if constexpr (std::is_convertible_v) { + return '"' + std::string{std::string_view{val}} + '"'; + } else if constexpr (IsContainer) { + std::stringstream out; + out << "{"; + const auto size = val.size(); + if (size > 1) { + out << repr(*val.begin()); + std::for_each( + std::next(val.begin()), + std::next( + val.begin(), + static_cast( + std::min(size, repr_max_container_size) - 1)), + [&out](const auto &v) { out << " " << repr(v); }); + if (size <= repr_max_container_size) { + out << " "; + } else { + out << "..."; + } + } + if (size > 0) { + out << repr(*std::prev(val.end())); + } + out << "}"; + return out.str(); + } else if constexpr (IsStreamable) { + std::stringstream out; + out << val; + return out.str(); + } else { + return ""; + } +} + +namespace { + +template constexpr bool standard_signed_integer = false; +template <> constexpr bool standard_signed_integer = true; +template <> constexpr bool standard_signed_integer = true; +template <> constexpr bool standard_signed_integer = true; +template <> constexpr bool standard_signed_integer = true; +template <> constexpr bool standard_signed_integer = true; + +template constexpr bool standard_unsigned_integer = false; +template <> constexpr bool standard_unsigned_integer = true; +template <> constexpr bool standard_unsigned_integer = true; +template <> constexpr bool standard_unsigned_integer = true; +template <> constexpr bool standard_unsigned_integer = true; +template <> +constexpr bool standard_unsigned_integer = true; + +} // namespace + +constexpr int radix_2 = 2; +constexpr int radix_8 = 8; +constexpr int radix_10 = 10; +constexpr int radix_16 = 16; + +template +constexpr bool standard_integer = + standard_signed_integer || standard_unsigned_integer; + +template +constexpr decltype(auto) +apply_plus_one_impl(F &&f, Tuple &&t, Extra &&x, + std::index_sequence /*unused*/) { + return std::invoke(std::forward(f), std::get(std::forward(t))..., + std::forward(x)); +} + +template +constexpr decltype(auto) apply_plus_one(F &&f, Tuple &&t, Extra &&x) { + return details::apply_plus_one_impl( + std::forward(f), std::forward(t), std::forward(x), + std::make_index_sequence< + std::tuple_size_v>>{}); +} + +constexpr auto pointer_range(std::string_view s) noexcept { + return std::tuple(s.data(), s.data() + s.size()); +} + +template +constexpr bool starts_with(std::basic_string_view prefix, + std::basic_string_view s) noexcept { + return s.substr(0, prefix.size()) == prefix; +} + +enum class chars_format { + scientific = 0xf1, + fixed = 0xf2, + hex = 0xf4, + binary = 0xf8, + general = fixed | scientific +}; + +struct ConsumeBinaryPrefixResult { + bool is_binary; + std::string_view rest; +}; + +constexpr auto consume_binary_prefix(std::string_view s) + -> ConsumeBinaryPrefixResult { + if (starts_with(std::string_view{"0b"}, s) || + starts_with(std::string_view{"0B"}, s)) { + s.remove_prefix(2); + return {true, s}; + } + return {false, s}; +} + +struct ConsumeHexPrefixResult { + bool is_hexadecimal; + std::string_view rest; +}; + +using namespace std::literals; + +constexpr auto consume_hex_prefix(std::string_view s) + -> ConsumeHexPrefixResult { + if (starts_with("0x"sv, s) || starts_with("0X"sv, s)) { + s.remove_prefix(2); + return {true, s}; + } + return {false, s}; +} + +template +inline auto do_from_chars(std::string_view s) -> T { + T x{0}; + auto [first, last] = pointer_range(s); + auto [ptr, ec] = std::from_chars(first, last, x, Param); + if (ec == std::errc()) { + if (ptr == last) { + return x; + } + throw std::invalid_argument{"pattern '" + std::string(s) + + "' does not match to the end"}; + } + if (ec == std::errc::invalid_argument) { + throw std::invalid_argument{"pattern '" + std::string(s) + "' not found"}; + } + if (ec == std::errc::result_out_of_range) { + throw std::range_error{"'" + std::string(s) + "' not representable"}; + } + return x; // unreachable +} + +template struct parse_number { + auto operator()(std::string_view s) -> T { + return do_from_chars(s); + } +}; + +template struct parse_number { + auto operator()(std::string_view s) -> T { + if (auto [ok, rest] = consume_binary_prefix(s); ok) { + return do_from_chars(rest); + } + throw std::invalid_argument{"pattern not found"}; + } +}; + +template struct parse_number { + auto operator()(std::string_view s) -> T { + if (starts_with("0x"sv, s) || starts_with("0X"sv, s)) { + if (auto [ok, rest] = consume_hex_prefix(s); ok) { + try { + return do_from_chars(rest); + } catch (const std::invalid_argument &err) { + throw std::invalid_argument("Failed to parse '" + std::string(s) + + "' as hexadecimal: " + err.what()); + } catch (const std::range_error &err) { + throw std::range_error("Failed to parse '" + std::string(s) + + "' as hexadecimal: " + err.what()); + } + } + } else { + // Allow passing hex numbers without prefix + // Shape 'x' already has to be specified + try { + return do_from_chars(s); + } catch (const std::invalid_argument &err) { + throw std::invalid_argument("Failed to parse '" + std::string(s) + + "' as hexadecimal: " + err.what()); + } catch (const std::range_error &err) { + throw std::range_error("Failed to parse '" + std::string(s) + + "' as hexadecimal: " + err.what()); + } + } + + throw std::invalid_argument{"pattern '" + std::string(s) + + "' not identified as hexadecimal"}; + } +}; + +template struct parse_number { + auto operator()(std::string_view s) -> T { + auto [ok, rest] = consume_hex_prefix(s); + if (ok) { + try { + return do_from_chars(rest); + } catch (const std::invalid_argument &err) { + throw std::invalid_argument("Failed to parse '" + std::string(s) + + "' as hexadecimal: " + err.what()); + } catch (const std::range_error &err) { + throw std::range_error("Failed to parse '" + std::string(s) + + "' as hexadecimal: " + err.what()); + } + } + + auto [ok_binary, rest_binary] = consume_binary_prefix(s); + if (ok_binary) { + try { + return do_from_chars(rest_binary); + } catch (const std::invalid_argument &err) { + throw std::invalid_argument("Failed to parse '" + std::string(s) + + "' as binary: " + err.what()); + } catch (const std::range_error &err) { + throw std::range_error("Failed to parse '" + std::string(s) + + "' as binary: " + err.what()); + } + } + + if (starts_with("0"sv, s)) { + try { + return do_from_chars(rest); + } catch (const std::invalid_argument &err) { + throw std::invalid_argument("Failed to parse '" + std::string(s) + + "' as octal: " + err.what()); + } catch (const std::range_error &err) { + throw std::range_error("Failed to parse '" + std::string(s) + + "' as octal: " + err.what()); + } + } + + try { + return do_from_chars(rest); + } catch (const std::invalid_argument &err) { + throw std::invalid_argument("Failed to parse '" + std::string(s) + + "' as decimal integer: " + err.what()); + } catch (const std::range_error &err) { + throw std::range_error("Failed to parse '" + std::string(s) + + "' as decimal integer: " + err.what()); + } + } +}; + +namespace { + +template inline const auto generic_strtod = nullptr; +template <> inline const auto generic_strtod = ARGPARSE_CUSTOM_STRTOF; +template <> inline const auto generic_strtod = ARGPARSE_CUSTOM_STRTOD; +template <> +inline const auto generic_strtod = ARGPARSE_CUSTOM_STRTOLD; + +} // namespace + +template inline auto do_strtod(std::string const &s) -> T { + if (isspace(static_cast(s[0])) || s[0] == '+') { + throw std::invalid_argument{"pattern '" + s + "' not found"}; + } + + auto [first, last] = pointer_range(s); + char *ptr; + + errno = 0; + auto x = generic_strtod(first, &ptr); + if (errno == 0) { + if (ptr == last) { + return x; + } + throw std::invalid_argument{"pattern '" + s + + "' does not match to the end"}; + } + if (errno == ERANGE) { + throw std::range_error{"'" + s + "' not representable"}; + } + return x; // unreachable +} + +template struct parse_number { + auto operator()(std::string const &s) -> T { + if (auto r = consume_hex_prefix(s); r.is_hexadecimal) { + throw std::invalid_argument{ + "chars_format::general does not parse hexfloat"}; + } + if (auto r = consume_binary_prefix(s); r.is_binary) { + throw std::invalid_argument{ + "chars_format::general does not parse binfloat"}; + } + + try { + return do_strtod(s); + } catch (const std::invalid_argument &err) { + throw std::invalid_argument("Failed to parse '" + s + + "' as number: " + err.what()); + } catch (const std::range_error &err) { + throw std::range_error("Failed to parse '" + s + + "' as number: " + err.what()); + } + } +}; + +template struct parse_number { + auto operator()(std::string const &s) -> T { + if (auto r = consume_hex_prefix(s); !r.is_hexadecimal) { + throw std::invalid_argument{"chars_format::hex parses hexfloat"}; + } + if (auto r = consume_binary_prefix(s); r.is_binary) { + throw std::invalid_argument{"chars_format::hex does not parse binfloat"}; + } + + try { + return do_strtod(s); + } catch (const std::invalid_argument &err) { + throw std::invalid_argument("Failed to parse '" + s + + "' as hexadecimal: " + err.what()); + } catch (const std::range_error &err) { + throw std::range_error("Failed to parse '" + s + + "' as hexadecimal: " + err.what()); + } + } +}; + +template struct parse_number { + auto operator()(std::string const &s) -> T { + if (auto r = consume_hex_prefix(s); r.is_hexadecimal) { + throw std::invalid_argument{ + "chars_format::binary does not parse hexfloat"}; + } + if (auto r = consume_binary_prefix(s); !r.is_binary) { + throw std::invalid_argument{"chars_format::binary parses binfloat"}; + } + + return do_strtod(s); + } +}; + +template struct parse_number { + auto operator()(std::string const &s) -> T { + if (auto r = consume_hex_prefix(s); r.is_hexadecimal) { + throw std::invalid_argument{ + "chars_format::scientific does not parse hexfloat"}; + } + if (auto r = consume_binary_prefix(s); r.is_binary) { + throw std::invalid_argument{ + "chars_format::scientific does not parse binfloat"}; + } + if (s.find_first_of("eE") == std::string::npos) { + throw std::invalid_argument{ + "chars_format::scientific requires exponent part"}; + } + + try { + return do_strtod(s); + } catch (const std::invalid_argument &err) { + throw std::invalid_argument("Failed to parse '" + s + + "' as scientific notation: " + err.what()); + } catch (const std::range_error &err) { + throw std::range_error("Failed to parse '" + s + + "' as scientific notation: " + err.what()); + } + } +}; + +template struct parse_number { + auto operator()(std::string const &s) -> T { + if (auto r = consume_hex_prefix(s); r.is_hexadecimal) { + throw std::invalid_argument{ + "chars_format::fixed does not parse hexfloat"}; + } + if (auto r = consume_binary_prefix(s); r.is_binary) { + throw std::invalid_argument{ + "chars_format::fixed does not parse binfloat"}; + } + if (s.find_first_of("eE") != std::string::npos) { + throw std::invalid_argument{ + "chars_format::fixed does not parse exponent part"}; + } + + try { + return do_strtod(s); + } catch (const std::invalid_argument &err) { + throw std::invalid_argument("Failed to parse '" + s + + "' as fixed notation: " + err.what()); + } catch (const std::range_error &err) { + throw std::range_error("Failed to parse '" + s + + "' as fixed notation: " + err.what()); + } + } +}; + +template +std::string join(StrIt first, StrIt last, const std::string &separator) { + if (first == last) { + return ""; + } + std::stringstream value; + value << *first; + ++first; + while (first != last) { + value << separator << *first; + ++first; + } + return value.str(); +} + +template struct can_invoke_to_string { + template + static auto test(int) + -> decltype(std::to_string(std::declval()), std::true_type{}); + + template static auto test(...) -> std::false_type; + + static constexpr bool value = decltype(test(0))::value; +}; + +template struct IsChoiceTypeSupported { + using CleanType = typename std::decay::type; + static const bool value = std::is_integral::value || + std::is_same::value || + std::is_same::value || + std::is_same::value; +}; + +template +std::size_t get_levenshtein_distance(const StringType &s1, + const StringType &s2) { + std::vector> dp( + s1.size() + 1, std::vector(s2.size() + 1, 0)); + + for (std::size_t i = 0; i <= s1.size(); ++i) { + for (std::size_t j = 0; j <= s2.size(); ++j) { + if (i == 0) { + dp[i][j] = j; + } else if (j == 0) { + dp[i][j] = i; + } else if (s1[i - 1] == s2[j - 1]) { + dp[i][j] = dp[i - 1][j - 1]; + } else { + dp[i][j] = 1 + std::min({dp[i - 1][j], dp[i][j - 1], dp[i - 1][j - 1]}); + } + } + } + + return dp[s1.size()][s2.size()]; +} + +template +std::string get_most_similar_string(const std::map &map, + const std::string &input) { + std::string most_similar{}; + std::size_t min_distance = (std::numeric_limits::max)(); + + for (const auto &entry : map) { + std::size_t distance = get_levenshtein_distance(entry.first, input); + if (distance < min_distance) { + min_distance = distance; + most_similar = entry.first; + } + } + + return most_similar; +} + +} // namespace details + +enum class nargs_pattern { optional, any, at_least_one }; + +enum class default_arguments : unsigned int { + none = 0, + help = 1, + version = 2, + all = help | version, +}; + +inline default_arguments operator&(const default_arguments &a, + const default_arguments &b) { + return static_cast( + static_cast::type>(a) & + static_cast::type>(b)); +} + +class ArgumentParser; + +class Argument { + friend class ArgumentParser; + friend auto operator<<(std::ostream &stream, const ArgumentParser &parser) + -> std::ostream &; + + template + explicit Argument(std::string_view prefix_chars, + std::array &&a, + std::index_sequence /*unused*/) + : m_accepts_optional_like_value(false), + m_is_optional((is_optional(a[I], prefix_chars) || ...)), + m_is_required(false), m_is_repeatable(false), m_is_used(false), + m_is_hidden(false), m_prefix_chars(prefix_chars) { + ((void)m_names.emplace_back(a[I]), ...); + std::sort( + m_names.begin(), m_names.end(), [](const auto &lhs, const auto &rhs) { + return lhs.size() == rhs.size() ? lhs < rhs : lhs.size() < rhs.size(); + }); + } + +public: + template + explicit Argument(std::string_view prefix_chars, + std::array &&a) + : Argument(prefix_chars, std::move(a), std::make_index_sequence{}) {} + + Argument &help(std::string help_text) { + m_help = std::move(help_text); + return *this; + } + + Argument &metavar(std::string metavar) { + m_metavar = std::move(metavar); + return *this; + } + + template Argument &default_value(T &&value) { + m_num_args_range = NArgsRange{0, m_num_args_range.get_max()}; + m_default_value_repr = details::repr(value); + + if constexpr (std::is_convertible_v) { + m_default_value_str = std::string{std::string_view{value}}; + } else if constexpr (details::can_invoke_to_string::value) { + m_default_value_str = std::to_string(value); + } + + m_default_value = std::forward(value); + return *this; + } + + Argument &default_value(const char *value) { + return default_value(std::string(value)); + } + + Argument &required() { + m_is_required = true; + return *this; + } + + Argument &implicit_value(std::any value) { + m_implicit_value = std::move(value); + m_num_args_range = NArgsRange{0, 0}; + return *this; + } + + // This is shorthand for: + // program.add_argument("foo") + // .default_value(false) + // .implicit_value(true) + Argument &flag() { + default_value(false); + implicit_value(true); + return *this; + } + + template + auto action(F &&callable, Args &&... bound_args) + -> std::enable_if_t, + Argument &> { + using action_type = std::conditional_t< + std::is_void_v>, + void_action, valued_action>; + if constexpr (sizeof...(Args) == 0) { + m_actions.emplace_back(std::forward(callable)); + } else { + m_actions.emplace_back( + [f = std::forward(callable), + tup = std::make_tuple(std::forward(bound_args)...)]( + std::string const &opt) mutable { + return details::apply_plus_one(f, tup, opt); + }); + } + return *this; + } + + auto &store_into(bool &var) { + if ((!m_default_value.has_value()) && (!m_implicit_value.has_value())) { + flag(); + } + if (m_default_value.has_value()) { + var = std::any_cast(m_default_value); + } + action([&var](const auto & /*unused*/) { + var = true; + return var; + }); + return *this; + } + + template ::value>::type * = nullptr> + auto &store_into(T &var) { + if (m_default_value.has_value()) { + var = std::any_cast(m_default_value); + } + action([&var](const auto &s) { + var = details::parse_number()(s); + return var; + }); + return *this; + } + + template ::value>::type * = nullptr> + auto &store_into(T &var) { + if (m_default_value.has_value()) { + var = std::any_cast(m_default_value); + } + action([&var](const auto &s) { + var = details::parse_number()(s); + return var; + }); + return *this; + } + + auto &store_into(std::string &var) { + if (m_default_value.has_value()) { + var = std::any_cast(m_default_value); + } + action([&var](const std::string &s) { + var = s; + return var; + }); + return *this; + } + + auto &store_into(std::filesystem::path &var) { + if (m_default_value.has_value()) { + var = std::any_cast(m_default_value); + } + action([&var](const std::string &s) { var = s; }); + return *this; + } + + auto &store_into(std::vector &var) { + if (m_default_value.has_value()) { + var = std::any_cast>(m_default_value); + } + action([this, &var](const std::string &s) { + if (!m_is_used) { + var.clear(); + } + m_is_used = true; + var.push_back(s); + return var; + }); + return *this; + } + + auto &store_into(std::vector &var) { + if (m_default_value.has_value()) { + var = std::any_cast>(m_default_value); + } + action([this, &var](const std::string &s) { + if (!m_is_used) { + var.clear(); + } + m_is_used = true; + var.push_back(details::parse_number()(s)); + return var; + }); + return *this; + } + + auto &store_into(std::set &var) { + if (m_default_value.has_value()) { + var = std::any_cast>(m_default_value); + } + action([this, &var](const std::string &s) { + if (!m_is_used) { + var.clear(); + } + m_is_used = true; + var.insert(s); + return var; + }); + return *this; + } + + auto &store_into(std::set &var) { + if (m_default_value.has_value()) { + var = std::any_cast>(m_default_value); + } + action([this, &var](const std::string &s) { + if (!m_is_used) { + var.clear(); + } + m_is_used = true; + var.insert(details::parse_number()(s)); + return var; + }); + return *this; + } + + auto &append() { + m_is_repeatable = true; + return *this; + } + + // Cause the argument to be invisible in usage and help + auto &hidden() { + m_is_hidden = true; + return *this; + } + + template + auto scan() -> std::enable_if_t, Argument &> { + static_assert(!(std::is_const_v || std::is_volatile_v), + "T should not be cv-qualified"); + auto is_one_of = [](char c, auto... x) constexpr { + return ((c == x) || ...); + }; + + if constexpr (is_one_of(Shape, 'd') && details::standard_integer) { + action(details::parse_number()); + } else if constexpr (is_one_of(Shape, 'i') && + details::standard_integer) { + action(details::parse_number()); + } else if constexpr (is_one_of(Shape, 'u') && + details::standard_unsigned_integer) { + action(details::parse_number()); + } else if constexpr (is_one_of(Shape, 'b') && + details::standard_unsigned_integer) { + action(details::parse_number()); + } else if constexpr (is_one_of(Shape, 'o') && + details::standard_unsigned_integer) { + action(details::parse_number()); + } else if constexpr (is_one_of(Shape, 'x', 'X') && + details::standard_unsigned_integer) { + action(details::parse_number()); + } else if constexpr (is_one_of(Shape, 'a', 'A') && + std::is_floating_point_v) { + action(details::parse_number()); + } else if constexpr (is_one_of(Shape, 'e', 'E') && + std::is_floating_point_v) { + action(details::parse_number()); + } else if constexpr (is_one_of(Shape, 'f', 'F') && + std::is_floating_point_v) { + action(details::parse_number()); + } else if constexpr (is_one_of(Shape, 'g', 'G') && + std::is_floating_point_v) { + action(details::parse_number()); + } else { + static_assert(alignof(T) == 0, "No scan specification for T"); + } + + return *this; + } + + Argument &nargs(std::size_t num_args) { + m_num_args_range = NArgsRange{num_args, num_args}; + return *this; + } + + Argument &nargs(std::size_t num_args_min, std::size_t num_args_max) { + m_num_args_range = NArgsRange{num_args_min, num_args_max}; + return *this; + } + + Argument &nargs(nargs_pattern pattern) { + switch (pattern) { + case nargs_pattern::optional: + m_num_args_range = NArgsRange{0, 1}; + break; + case nargs_pattern::any: + m_num_args_range = + NArgsRange{0, (std::numeric_limits::max)()}; + break; + case nargs_pattern::at_least_one: + m_num_args_range = + NArgsRange{1, (std::numeric_limits::max)()}; + break; + } + return *this; + } + + Argument &remaining() { + m_accepts_optional_like_value = true; + return nargs(nargs_pattern::any); + } + + template void add_choice(T &&choice) { + static_assert(details::IsChoiceTypeSupported::value, + "Only string or integer type supported for choice"); + static_assert(std::is_convertible_v || + details::can_invoke_to_string::value, + "Choice is not convertible to string_type"); + if (!m_choices.has_value()) { + m_choices = std::vector{}; + } + + if constexpr (std::is_convertible_v) { + m_choices.value().push_back( + std::string{std::string_view{std::forward(choice)}}); + } else if constexpr (details::can_invoke_to_string::value) { + m_choices.value().push_back(std::to_string(std::forward(choice))); + } + } + + Argument &choices() { + if (!m_choices.has_value()) { + throw std::runtime_error("Zero choices provided"); + } + return *this; + } + + template + Argument &choices(T &&first, U &&... rest) { + add_choice(std::forward(first)); + choices(std::forward(rest)...); + return *this; + } + + void find_default_value_in_choices_or_throw() const { + + const auto &choices = m_choices.value(); + + if (m_default_value.has_value()) { + if (std::find(choices.begin(), choices.end(), m_default_value_str) == + choices.end()) { + // provided arg not in list of allowed choices + // report error + + std::string choices_as_csv = + std::accumulate(choices.begin(), choices.end(), std::string(), + [](const std::string &a, const std::string &b) { + return a + (a.empty() ? "" : ", ") + b; + }); + + throw std::runtime_error( + std::string{"Invalid default value "} + m_default_value_repr + + " - allowed options: {" + choices_as_csv + "}"); + } + } + } + + template + bool is_value_in_choices(Iterator option_it) const { + + const auto &choices = m_choices.value(); + + return (std::find(choices.begin(), choices.end(), *option_it) != + choices.end()); + } + + template + void throw_invalid_arguments_error(Iterator option_it) const { + const auto &choices = m_choices.value(); + const std::string choices_as_csv = std::accumulate( + choices.begin(), choices.end(), std::string(), + [](const std::string &option_a, const std::string &option_b) { + return option_a + (option_a.empty() ? "" : ", ") + option_b; + }); + + throw std::runtime_error(std::string{"Invalid argument "} + + details::repr(*option_it) + + " - allowed options: {" + choices_as_csv + "}"); + } + + /* The dry_run parameter can be set to true to avoid running the actions, + * and setting m_is_used. This may be used by a pre-processing step to do + * a first iteration over arguments. + */ + template + Iterator consume(Iterator start, Iterator end, + std::string_view used_name = {}, bool dry_run = false) { + if (!m_is_repeatable && m_is_used) { + throw std::runtime_error( + std::string("Duplicate argument ").append(used_name)); + } + m_used_name = used_name; + + std::size_t passed_options = 0; + + if (m_choices.has_value()) { + // Check each value in (start, end) and make sure + // it is in the list of allowed choices/options + const auto max_number_of_args = m_num_args_range.get_max(); + const auto min_number_of_args = m_num_args_range.get_min(); + for (auto it = start; it != end; ++it) { + if (is_value_in_choices(it)) { + passed_options += 1; + continue; + } + + if ((passed_options >= min_number_of_args) && + (passed_options <= max_number_of_args)) { + break; + } + + throw_invalid_arguments_error(it); + } + } + + const auto num_args_max = + (m_choices.has_value()) ? passed_options : m_num_args_range.get_max(); + const auto num_args_min = m_num_args_range.get_min(); + std::size_t dist = 0; + if (num_args_max == 0) { + if (!dry_run) { + m_values.emplace_back(m_implicit_value); + for(auto &action: m_actions) { + std::visit([&](const auto &f) { f({}); }, action); + } + if(m_actions.empty()){ + std::visit([&](const auto &f) { f({}); }, m_default_action); + } + m_is_used = true; + } + return start; + } + if ((dist = static_cast(std::distance(start, end))) >= + num_args_min) { + if (num_args_max < dist) { + end = std::next(start, static_cast( + num_args_max)); + } + if (!m_accepts_optional_like_value) { + end = std::find_if( + start, end, + std::bind(is_optional, std::placeholders::_1, m_prefix_chars)); + dist = static_cast(std::distance(start, end)); + if (dist < num_args_min) { + throw std::runtime_error("Too few arguments for '" + + std::string(m_used_name) + "'."); + } + } + struct ActionApply { + void operator()(valued_action &f) { + std::transform(first, last, std::back_inserter(self.m_values), f); + } + + void operator()(void_action &f) { + std::for_each(first, last, f); + if (!self.m_default_value.has_value()) { + if (!self.m_accepts_optional_like_value) { + self.m_values.resize( + static_cast(std::distance(first, last))); + } + } + } + + Iterator first, last; + Argument &self; + }; + if (!dry_run) { + for(auto &action: m_actions) { + std::visit(ActionApply{start, end, *this}, action); + } + if(m_actions.empty()){ + std::visit(ActionApply{start, end, *this}, m_default_action); + } + m_is_used = true; + } + return end; + } + if (m_default_value.has_value()) { + if (!dry_run) { + m_is_used = true; + } + return start; + } + throw std::runtime_error("Too few arguments for '" + + std::string(m_used_name) + "'."); + } + + /* + * @throws std::runtime_error if argument values are not valid + */ + void validate() const { + if (m_is_optional) { + // TODO: check if an implicit value was programmed for this argument + if (!m_is_used && !m_default_value.has_value() && m_is_required) { + throw_required_arg_not_used_error(); + } + if (m_is_used && m_is_required && m_values.empty()) { + throw_required_arg_no_value_provided_error(); + } + } else { + if (!m_num_args_range.contains(m_values.size()) && + !m_default_value.has_value()) { + throw_nargs_range_validation_error(); + } + } + + if (m_choices.has_value()) { + // Make sure the default value (if provided) + // is in the list of choices + find_default_value_in_choices_or_throw(); + } + } + + std::string get_names_csv(char separator = ',') const { + return std::accumulate( + m_names.begin(), m_names.end(), std::string{""}, + [&](const std::string &result, const std::string &name) { + return result.empty() ? name : result + separator + name; + }); + } + + std::string get_usage_full() const { + std::stringstream usage; + + usage << get_names_csv('/'); + const std::string metavar = !m_metavar.empty() ? m_metavar : "VAR"; + if (m_num_args_range.get_max() > 0) { + usage << " " << metavar; + if (m_num_args_range.get_max() > 1) { + usage << "..."; + } + } + return usage.str(); + } + + std::string get_inline_usage() const { + std::stringstream usage; + // Find the longest variant to show in the usage string + std::string longest_name = m_names.front(); + for (const auto &s : m_names) { + if (s.size() > longest_name.size()) { + longest_name = s; + } + } + if (!m_is_required) { + usage << "["; + } + usage << longest_name; + const std::string metavar = !m_metavar.empty() ? m_metavar : "VAR"; + if (m_num_args_range.get_max() > 0) { + usage << " " << metavar; + if (m_num_args_range.get_max() > 1 && + m_metavar.find("> <") == std::string::npos) { + usage << "..."; + } + } + if (!m_is_required) { + usage << "]"; + } + if (m_is_repeatable) { + usage << "..."; + } + return usage.str(); + } + + std::size_t get_arguments_length() const { + + std::size_t names_size = std::accumulate( + std::begin(m_names), std::end(m_names), std::size_t(0), + [](const auto &sum, const auto &s) { return sum + s.size(); }); + + if (is_positional(m_names.front(), m_prefix_chars)) { + // A set metavar means this replaces the names + if (!m_metavar.empty()) { + // Indent and metavar + return 2 + m_metavar.size(); + } + + // Indent and space-separated + return 2 + names_size + (m_names.size() - 1); + } + // Is an option - include both names _and_ metavar + // size = text + (", " between names) + std::size_t size = names_size + 2 * (m_names.size() - 1); + if (!m_metavar.empty() && m_num_args_range == NArgsRange{1, 1}) { + size += m_metavar.size() + 1; + } + return size + 2; // indent + } + + friend std::ostream &operator<<(std::ostream &stream, + const Argument &argument) { + std::stringstream name_stream; + name_stream << " "; // indent + if (argument.is_positional(argument.m_names.front(), + argument.m_prefix_chars)) { + if (!argument.m_metavar.empty()) { + name_stream << argument.m_metavar; + } else { + name_stream << details::join(argument.m_names.begin(), + argument.m_names.end(), " "); + } + } else { + name_stream << details::join(argument.m_names.begin(), + argument.m_names.end(), ", "); + // If we have a metavar, and one narg - print the metavar + if (!argument.m_metavar.empty() && + argument.m_num_args_range == NArgsRange{1, 1}) { + name_stream << " " << argument.m_metavar; + } + else if (!argument.m_metavar.empty() && + argument.m_num_args_range.get_min() == argument.m_num_args_range.get_max() && + argument.m_metavar.find("> <") != std::string::npos) { + name_stream << " " << argument.m_metavar; + } + } + + // align multiline help message + auto stream_width = stream.width(); + auto name_padding = std::string(name_stream.str().size(), ' '); + auto pos = std::string::size_type{}; + auto prev = std::string::size_type{}; + auto first_line = true; + auto hspace = " "; // minimal space between name and help message + stream << name_stream.str(); + std::string_view help_view(argument.m_help); + while ((pos = argument.m_help.find('\n', prev)) != std::string::npos) { + auto line = help_view.substr(prev, pos - prev + 1); + if (first_line) { + stream << hspace << line; + first_line = false; + } else { + stream.width(stream_width); + stream << name_padding << hspace << line; + } + prev += pos - prev + 1; + } + if (first_line) { + stream << hspace << argument.m_help; + } else { + auto leftover = help_view.substr(prev, argument.m_help.size() - prev); + if (!leftover.empty()) { + stream.width(stream_width); + stream << name_padding << hspace << leftover; + } + } + + // print nargs spec + if (!argument.m_help.empty()) { + stream << " "; + } + stream << argument.m_num_args_range; + + bool add_space = false; + if (argument.m_default_value.has_value() && + argument.m_num_args_range != NArgsRange{0, 0}) { + stream << "[default: " << argument.m_default_value_repr << "]"; + add_space = true; + } else if (argument.m_is_required) { + stream << "[required]"; + add_space = true; + } + if (argument.m_is_repeatable) { + if (add_space) { + stream << " "; + } + stream << "[may be repeated]"; + } + stream << "\n"; + return stream; + } + + template bool operator!=(const T &rhs) const { + return !(*this == rhs); + } + + /* + * Compare to an argument value of known type + * @throws std::logic_error in case of incompatible types + */ + template bool operator==(const T &rhs) const { + if constexpr (!details::IsContainer) { + return get() == rhs; + } else { + using ValueType = typename T::value_type; + auto lhs = get(); + return std::equal(std::begin(lhs), std::end(lhs), std::begin(rhs), + std::end(rhs), [](const auto &a, const auto &b) { + return std::any_cast(a) == b; + }); + } + } + + /* + * positional: + * _empty_ + * '-' + * '-' decimal-literal + * !'-' anything + */ + static bool is_positional(std::string_view name, + std::string_view prefix_chars) { + auto first = lookahead(name); + + if (first == eof) { + return true; + } + if (prefix_chars.find(static_cast(first)) != + std::string_view::npos) { + name.remove_prefix(1); + if (name.empty()) { + return true; + } + return is_decimal_literal(name); + } + return true; + } + +private: + class NArgsRange { + std::size_t m_min; + std::size_t m_max; + + public: + NArgsRange(std::size_t minimum, std::size_t maximum) + : m_min(minimum), m_max(maximum) { + if (minimum > maximum) { + throw std::logic_error("Range of number of arguments is invalid"); + } + } + + bool contains(std::size_t value) const { + return value >= m_min && value <= m_max; + } + + bool is_exact() const { return m_min == m_max; } + + bool is_right_bounded() const { + return m_max < (std::numeric_limits::max)(); + } + + std::size_t get_min() const { return m_min; } + + std::size_t get_max() const { return m_max; } + + // Print help message + friend auto operator<<(std::ostream &stream, const NArgsRange &range) + -> std::ostream & { + if (range.m_min == range.m_max) { + if (range.m_min != 0 && range.m_min != 1) { + stream << "[nargs: " << range.m_min << "] "; + } + } else { + if (range.m_max == (std::numeric_limits::max)()) { + stream << "[nargs: " << range.m_min << " or more] "; + } else { + stream << "[nargs=" << range.m_min << ".." << range.m_max << "] "; + } + } + return stream; + } + + bool operator==(const NArgsRange &rhs) const { + return rhs.m_min == m_min && rhs.m_max == m_max; + } + + bool operator!=(const NArgsRange &rhs) const { return !(*this == rhs); } + }; + + void throw_nargs_range_validation_error() const { + std::stringstream stream; + if (!m_used_name.empty()) { + stream << m_used_name << ": "; + } else { + stream << m_names.front() << ": "; + } + if (m_num_args_range.is_exact()) { + stream << m_num_args_range.get_min(); + } else if (m_num_args_range.is_right_bounded()) { + stream << m_num_args_range.get_min() << " to " + << m_num_args_range.get_max(); + } else { + stream << m_num_args_range.get_min() << " or more"; + } + stream << " argument(s) expected. " << m_values.size() << " provided."; + throw std::runtime_error(stream.str()); + } + + void throw_required_arg_not_used_error() const { + std::stringstream stream; + stream << m_names.front() << ": required."; + throw std::runtime_error(stream.str()); + } + + void throw_required_arg_no_value_provided_error() const { + std::stringstream stream; + stream << m_used_name << ": no value provided."; + throw std::runtime_error(stream.str()); + } + + static constexpr int eof = std::char_traits::eof(); + + static auto lookahead(std::string_view s) -> int { + if (s.empty()) { + return eof; + } + return static_cast(static_cast(s[0])); + } + + /* + * decimal-literal: + * '0' + * nonzero-digit digit-sequence_opt + * integer-part fractional-part + * fractional-part + * integer-part '.' exponent-part_opt + * integer-part exponent-part + * + * integer-part: + * digit-sequence + * + * fractional-part: + * '.' post-decimal-point + * + * post-decimal-point: + * digit-sequence exponent-part_opt + * + * exponent-part: + * 'e' post-e + * 'E' post-e + * + * post-e: + * sign_opt digit-sequence + * + * sign: one of + * '+' '-' + */ + static bool is_decimal_literal(std::string_view s) { + auto is_digit = [](auto c) constexpr { + switch (c) { + case '0': + case '1': + case '2': + case '3': + case '4': + case '5': + case '6': + case '7': + case '8': + case '9': + return true; + default: + return false; + } + }; + + // precondition: we have consumed or will consume at least one digit + auto consume_digits = [=](std::string_view sd) { + // NOLINTNEXTLINE(readability-qualified-auto) + auto it = std::find_if_not(std::begin(sd), std::end(sd), is_digit); + return sd.substr(static_cast(it - std::begin(sd))); + }; + + switch (lookahead(s)) { + case '0': { + s.remove_prefix(1); + if (s.empty()) { + return true; + } + goto integer_part; + } + case '1': + case '2': + case '3': + case '4': + case '5': + case '6': + case '7': + case '8': + case '9': { + s = consume_digits(s); + if (s.empty()) { + return true; + } + goto integer_part_consumed; + } + case '.': { + s.remove_prefix(1); + goto post_decimal_point; + } + default: + return false; + } + + integer_part: + s = consume_digits(s); + integer_part_consumed: + switch (lookahead(s)) { + case '.': { + s.remove_prefix(1); + if (is_digit(lookahead(s))) { + goto post_decimal_point; + } else { + goto exponent_part_opt; + } + } + case 'e': + case 'E': { + s.remove_prefix(1); + goto post_e; + } + default: + return false; + } + + post_decimal_point: + if (is_digit(lookahead(s))) { + s = consume_digits(s); + goto exponent_part_opt; + } + return false; + + exponent_part_opt: + switch (lookahead(s)) { + case eof: + return true; + case 'e': + case 'E': { + s.remove_prefix(1); + goto post_e; + } + default: + return false; + } + + post_e: + switch (lookahead(s)) { + case '-': + case '+': + s.remove_prefix(1); + } + if (is_digit(lookahead(s))) { + s = consume_digits(s); + return s.empty(); + } + return false; + } + + static bool is_optional(std::string_view name, + std::string_view prefix_chars) { + return !is_positional(name, prefix_chars); + } + + /* + * Get argument value given a type + * @throws std::logic_error in case of incompatible types + */ + template T get() const { + if (!m_values.empty()) { + if constexpr (details::IsContainer) { + return any_cast_container(m_values); + } else { + return std::any_cast(m_values.front()); + } + } + if (m_default_value.has_value()) { + return std::any_cast(m_default_value); + } + if constexpr (details::IsContainer) { + if (!m_accepts_optional_like_value) { + return any_cast_container(m_values); + } + } + + throw std::logic_error("No value provided for '" + m_names.back() + "'."); + } + + /* + * Get argument value given a type. + * @pre The object has no default value. + * @returns The stored value if any, std::nullopt otherwise. + */ + template auto present() const -> std::optional { + if (m_default_value.has_value()) { + throw std::logic_error("Argument with default value always presents"); + } + if (m_values.empty()) { + return std::nullopt; + } + if constexpr (details::IsContainer) { + return any_cast_container(m_values); + } + return std::any_cast(m_values.front()); + } + + template + static auto any_cast_container(const std::vector &operand) -> T { + using ValueType = typename T::value_type; + + T result; + std::transform( + std::begin(operand), std::end(operand), std::back_inserter(result), + [](const auto &value) { return std::any_cast(value); }); + return result; + } + + void set_usage_newline_counter(int i) { m_usage_newline_counter = i; } + + void set_group_idx(std::size_t i) { m_group_idx = i; } + + std::vector m_names; + std::string_view m_used_name; + std::string m_help; + std::string m_metavar; + std::any m_default_value; + std::string m_default_value_repr; + std::optional + m_default_value_str; // used for checking default_value against choices + std::any m_implicit_value; + std::optional> m_choices{std::nullopt}; + using valued_action = std::function; + using void_action = std::function; + std::vector> m_actions; + std::variant m_default_action{ + std::in_place_type, + [](const std::string &value) { return value; }}; + std::vector m_values; + NArgsRange m_num_args_range{1, 1}; + // Bit field of bool values. Set default value in ctor. + bool m_accepts_optional_like_value : 1; + bool m_is_optional : 1; + bool m_is_required : 1; + bool m_is_repeatable : 1; + bool m_is_used : 1; + bool m_is_hidden : 1; // if set, does not appear in usage or help + std::string_view m_prefix_chars; // ArgumentParser has the prefix_chars + int m_usage_newline_counter = 0; + std::size_t m_group_idx = 0; +}; + +class ArgumentParser { +public: + explicit ArgumentParser(std::string program_name = {}, + std::string version = "1.0", + default_arguments add_args = default_arguments::all, + bool exit_on_default_arguments = true, + std::ostream &os = std::cout) + : m_program_name(std::move(program_name)), m_version(std::move(version)), + m_exit_on_default_arguments(exit_on_default_arguments), + m_parser_path(m_program_name) { + if ((add_args & default_arguments::help) == default_arguments::help) { + add_argument("-h", "--help") + .action([&](const auto & /*unused*/) { + os << help().str(); + if (m_exit_on_default_arguments) { + std::exit(0); + } + }) + .default_value(false) + .help("shows help message and exits") + .implicit_value(true) + .nargs(0); + } + if ((add_args & default_arguments::version) == default_arguments::version) { + add_argument("-v", "--version") + .action([&](const auto & /*unused*/) { + os << m_version << std::endl; + if (m_exit_on_default_arguments) { + std::exit(0); + } + }) + .default_value(false) + .help("prints version information and exits") + .implicit_value(true) + .nargs(0); + } + } + + ~ArgumentParser() = default; + + // ArgumentParser is meant to be used in a single function. + // Setup everything and parse arguments in one place. + // + // ArgumentParser internally uses std::string_views, + // references, iterators, etc. + // Many of these elements become invalidated after a copy or move. + ArgumentParser(const ArgumentParser &other) = delete; + ArgumentParser &operator=(const ArgumentParser &other) = delete; + ArgumentParser(ArgumentParser &&) noexcept = delete; + ArgumentParser &operator=(ArgumentParser &&) = delete; + + explicit operator bool() const { + auto arg_used = std::any_of(m_argument_map.cbegin(), m_argument_map.cend(), + [](auto &it) { return it.second->m_is_used; }); + auto subparser_used = + std::any_of(m_subparser_used.cbegin(), m_subparser_used.cend(), + [](auto &it) { return it.second; }); + + return m_is_parsed && (arg_used || subparser_used); + } + + // Parameter packing + // Call add_argument with variadic number of string arguments + template Argument &add_argument(Targs... f_args) { + using array_of_sv = std::array; + auto argument = + m_optional_arguments.emplace(std::cend(m_optional_arguments), + m_prefix_chars, array_of_sv{f_args...}); + + if (!argument->m_is_optional) { + m_positional_arguments.splice(std::cend(m_positional_arguments), + m_optional_arguments, argument); + } + argument->set_usage_newline_counter(m_usage_newline_counter); + argument->set_group_idx(m_group_names.size()); + + index_argument(argument); + return *argument; + } + + class MutuallyExclusiveGroup { + friend class ArgumentParser; + + public: + MutuallyExclusiveGroup() = delete; + + explicit MutuallyExclusiveGroup(ArgumentParser &parent, + bool required = false) + : m_parent(parent), m_required(required), m_elements({}) {} + + MutuallyExclusiveGroup(const MutuallyExclusiveGroup &other) = delete; + MutuallyExclusiveGroup & + operator=(const MutuallyExclusiveGroup &other) = delete; + + MutuallyExclusiveGroup(MutuallyExclusiveGroup &&other) noexcept + : m_parent(other.m_parent), m_required(other.m_required), + m_elements(std::move(other.m_elements)) { + other.m_elements.clear(); + } + + template Argument &add_argument(Targs... f_args) { + auto &argument = m_parent.add_argument(std::forward(f_args)...); + m_elements.push_back(&argument); + argument.set_usage_newline_counter(m_parent.m_usage_newline_counter); + argument.set_group_idx(m_parent.m_group_names.size()); + return argument; + } + + private: + ArgumentParser &m_parent; + bool m_required{false}; + std::vector m_elements{}; + }; + + MutuallyExclusiveGroup &add_mutually_exclusive_group(bool required = false) { + m_mutually_exclusive_groups.emplace_back(*this, required); + return m_mutually_exclusive_groups.back(); + } + + // Parameter packed add_parents method + // Accepts a variadic number of ArgumentParser objects + template + ArgumentParser &add_parents(const Targs &... f_args) { + for (const ArgumentParser &parent_parser : {std::ref(f_args)...}) { + for (const auto &argument : parent_parser.m_positional_arguments) { + auto it = m_positional_arguments.insert( + std::cend(m_positional_arguments), argument); + index_argument(it); + } + for (const auto &argument : parent_parser.m_optional_arguments) { + auto it = m_optional_arguments.insert(std::cend(m_optional_arguments), + argument); + index_argument(it); + } + } + return *this; + } + + // Ask for the next optional arguments to be displayed on a separate + // line in usage() output. Only effective if set_usage_max_line_width() is + // also used. + ArgumentParser &add_usage_newline() { + ++m_usage_newline_counter; + return *this; + } + + // Ask for the next optional arguments to be displayed in a separate section + // in usage() and help (<< *this) output. + // For usage(), this is only effective if set_usage_max_line_width() is + // also used. + ArgumentParser &add_group(std::string group_name) { + m_group_names.emplace_back(std::move(group_name)); + return *this; + } + + ArgumentParser &add_description(std::string description) { + m_description = std::move(description); + return *this; + } + + ArgumentParser &add_epilog(std::string epilog) { + m_epilog = std::move(epilog); + return *this; + } + + // Add a un-documented/hidden alias for an argument. + // Ideally we'd want this to be a method of Argument, but Argument + // does not own its owing ArgumentParser. + ArgumentParser &add_hidden_alias_for(Argument &arg, std::string_view alias) { + for (auto it = m_optional_arguments.begin(); + it != m_optional_arguments.end(); ++it) { + if (&(*it) == &arg) { + m_argument_map.insert_or_assign(std::string(alias), it); + return *this; + } + } + throw std::logic_error( + "Argument is not an optional argument of this parser"); + } + + /* Getter for arguments and subparsers. + * @throws std::logic_error in case of an invalid argument or subparser name + */ + template T &at(std::string_view name) { + if constexpr (std::is_same_v) { + return (*this)[name]; + } else { + std::string str_name(name); + auto subparser_it = m_subparser_map.find(str_name); + if (subparser_it != m_subparser_map.end()) { + return subparser_it->second->get(); + } + throw std::logic_error("No such subparser: " + str_name); + } + } + + ArgumentParser &set_prefix_chars(std::string prefix_chars) { + m_prefix_chars = std::move(prefix_chars); + return *this; + } + + ArgumentParser &set_assign_chars(std::string assign_chars) { + m_assign_chars = std::move(assign_chars); + return *this; + } + + /* Call parse_args_internal - which does all the work + * Then, validate the parsed arguments + * This variant is used mainly for testing + * @throws std::runtime_error in case of any invalid argument + */ + void parse_args(const std::vector &arguments) { + parse_args_internal(arguments); + // Check if all arguments are parsed + for ([[maybe_unused]] const auto &[unused, argument] : m_argument_map) { + argument->validate(); + } + + // Check each mutually exclusive group and make sure + // there are no constraint violations + for (const auto &group : m_mutually_exclusive_groups) { + auto mutex_argument_used{false}; + Argument *mutex_argument_it{nullptr}; + for (Argument *arg : group.m_elements) { + if (!mutex_argument_used && arg->m_is_used) { + mutex_argument_used = true; + mutex_argument_it = arg; + } else if (mutex_argument_used && arg->m_is_used) { + // Violation + throw std::runtime_error("Argument '" + arg->get_usage_full() + + "' not allowed with '" + + mutex_argument_it->get_usage_full() + "'"); + } + } + + if (!mutex_argument_used && group.m_required) { + // at least one argument from the group is + // required + std::string argument_names{}; + std::size_t i = 0; + std::size_t size = group.m_elements.size(); + for (Argument *arg : group.m_elements) { + if (i + 1 == size) { + // last + argument_names += std::string("'") + arg->get_usage_full() + std::string("' "); + } else { + argument_names += std::string("'") + arg->get_usage_full() + std::string("' or "); + } + i += 1; + } + throw std::runtime_error("One of the arguments " + argument_names + + "is required"); + } + } + } + + /* Call parse_known_args_internal - which does all the work + * Then, validate the parsed arguments + * This variant is used mainly for testing + * @throws std::runtime_error in case of any invalid argument + */ + std::vector + parse_known_args(const std::vector &arguments) { + auto unknown_arguments = parse_known_args_internal(arguments); + // Check if all arguments are parsed + for ([[maybe_unused]] const auto &[unused, argument] : m_argument_map) { + argument->validate(); + } + return unknown_arguments; + } + + /* Main entry point for parsing command-line arguments using this + * ArgumentParser + * @throws std::runtime_error in case of any invalid argument + */ + // NOLINTNEXTLINE(cppcoreguidelines-avoid-c-arrays) + void parse_args(int argc, const char *const argv[]) { + parse_args({argv, argv + argc}); + } + + /* Main entry point for parsing command-line arguments using this + * ArgumentParser + * @throws std::runtime_error in case of any invalid argument + */ + // NOLINTNEXTLINE(cppcoreguidelines-avoid-c-arrays) + auto parse_known_args(int argc, const char *const argv[]) { + return parse_known_args({argv, argv + argc}); + } + + /* Getter for options with default values. + * @throws std::logic_error if parse_args() has not been previously called + * @throws std::logic_error if there is no such option + * @throws std::logic_error if the option has no value + * @throws std::bad_any_cast if the option is not of type T + */ + template T get(std::string_view arg_name) const { + if (!m_is_parsed) { + throw std::logic_error("Nothing parsed, no arguments are available."); + } + return (*this)[arg_name].get(); + } + + /* Getter for options without default values. + * @pre The option has no default value. + * @throws std::logic_error if there is no such option + * @throws std::bad_any_cast if the option is not of type T + */ + template + auto present(std::string_view arg_name) const -> std::optional { + return (*this)[arg_name].present(); + } + + /* Getter that returns true for user-supplied options. Returns false if not + * user-supplied, even with a default value. + */ + auto is_used(std::string_view arg_name) const { + return (*this)[arg_name].m_is_used; + } + + /* Getter that returns true if a subcommand is used. + */ + auto is_subcommand_used(std::string_view subcommand_name) const { + return m_subparser_used.at(std::string(subcommand_name)); + } + + /* Getter that returns true if a subcommand is used. + */ + auto is_subcommand_used(const ArgumentParser &subparser) const { + return is_subcommand_used(subparser.m_program_name); + } + + /* Indexing operator. Return a reference to an Argument object + * Used in conjunction with Argument.operator== e.g., parser["foo"] == true + * @throws std::logic_error in case of an invalid argument name + */ + Argument &operator[](std::string_view arg_name) const { + std::string name(arg_name); + auto it = m_argument_map.find(name); + if (it != m_argument_map.end()) { + return *(it->second); + } + if (!is_valid_prefix_char(arg_name.front())) { + const auto legal_prefix_char = get_any_valid_prefix_char(); + const auto prefix = std::string(1, legal_prefix_char); + + // "-" + arg_name + name = prefix + name; + it = m_argument_map.find(name); + if (it != m_argument_map.end()) { + return *(it->second); + } + // "--" + arg_name + name = prefix + name; + it = m_argument_map.find(name); + if (it != m_argument_map.end()) { + return *(it->second); + } + } + throw std::logic_error("No such argument: " + std::string(arg_name)); + } + + // Print help message + friend auto operator<<(std::ostream &stream, const ArgumentParser &parser) + -> std::ostream & { + stream.setf(std::ios_base::left); + + auto longest_arg_length = parser.get_length_of_longest_argument(); + + stream << parser.usage() << "\n\n"; + + if (!parser.m_description.empty()) { + stream << parser.m_description << "\n\n"; + } + + const bool has_visible_positional_args = std::find_if( + parser.m_positional_arguments.begin(), + parser.m_positional_arguments.end(), + [](const auto &argument) { + return !argument.m_is_hidden; }) != + parser.m_positional_arguments.end(); + if (has_visible_positional_args) { + stream << "Positional arguments:\n"; + } + + for (const auto &argument : parser.m_positional_arguments) { + if (!argument.m_is_hidden) { + stream.width(static_cast(longest_arg_length)); + stream << argument; + } + } + + if (!parser.m_optional_arguments.empty()) { + stream << (!has_visible_positional_args ? "" : "\n") + << "Optional arguments:\n"; + } + + for (const auto &argument : parser.m_optional_arguments) { + if (argument.m_group_idx == 0 && !argument.m_is_hidden) { + stream.width(static_cast(longest_arg_length)); + stream << argument; + } + } + + for (size_t i_group = 0; i_group < parser.m_group_names.size(); ++i_group) { + stream << "\n" << parser.m_group_names[i_group] << " (detailed usage):\n"; + for (const auto &argument : parser.m_optional_arguments) { + if (argument.m_group_idx == i_group + 1 && !argument.m_is_hidden) { + stream.width(static_cast(longest_arg_length)); + stream << argument; + } + } + } + + bool has_visible_subcommands = std::any_of( + parser.m_subparser_map.begin(), parser.m_subparser_map.end(), + [](auto &p) { return !p.second->get().m_suppress; }); + + if (has_visible_subcommands) { + stream << (parser.m_positional_arguments.empty() + ? (parser.m_optional_arguments.empty() ? "" : "\n") + : "\n") + << "Subcommands:\n"; + for (const auto &[command, subparser] : parser.m_subparser_map) { + if (subparser->get().m_suppress) { + continue; + } + + stream << std::setw(2) << " "; + stream << std::setw(static_cast(longest_arg_length - 2)) + << command; + stream << " " << subparser->get().m_description << "\n"; + } + } + + if (!parser.m_epilog.empty()) { + stream << '\n'; + stream << parser.m_epilog << "\n\n"; + } + + return stream; + } + + // Format help message + auto help() const -> std::stringstream { + std::stringstream out; + out << *this; + return out; + } + + // Sets the maximum width for a line of the Usage message + ArgumentParser &set_usage_max_line_width(size_t w) { + this->m_usage_max_line_width = w; + return *this; + } + + // Asks to display arguments of mutually exclusive group on separate lines in + // the Usage message + ArgumentParser &set_usage_break_on_mutex() { + this->m_usage_break_on_mutex = true; + return *this; + } + + // Format usage part of help only + auto usage() const -> std::string { + std::stringstream stream; + + std::string curline("Usage: "); + curline += this->m_parser_path; + const bool multiline_usage = + this->m_usage_max_line_width < (std::numeric_limits::max)(); + const size_t indent_size = curline.size(); + + const auto deal_with_options_of_group = [&](std::size_t group_idx) { + bool found_options = false; + // Add any options inline here + const MutuallyExclusiveGroup *cur_mutex = nullptr; + int usage_newline_counter = -1; + for (const auto &argument : this->m_optional_arguments) { + if (argument.m_is_hidden) { + continue; + } + if (multiline_usage) { + if (argument.m_group_idx != group_idx) { + continue; + } + if (usage_newline_counter != argument.m_usage_newline_counter) { + if (usage_newline_counter >= 0) { + if (curline.size() > indent_size) { + stream << curline << std::endl; + curline = std::string(indent_size, ' '); + } + } + usage_newline_counter = argument.m_usage_newline_counter; + } + } + found_options = true; + const std::string arg_inline_usage = argument.get_inline_usage(); + const MutuallyExclusiveGroup *arg_mutex = + get_belonging_mutex(&argument); + if ((cur_mutex != nullptr) && (arg_mutex == nullptr)) { + curline += ']'; + if (this->m_usage_break_on_mutex) { + stream << curline << std::endl; + curline = std::string(indent_size, ' '); + } + } else if ((cur_mutex == nullptr) && (arg_mutex != nullptr)) { + if ((this->m_usage_break_on_mutex && curline.size() > indent_size) || + curline.size() + 3 + arg_inline_usage.size() > + this->m_usage_max_line_width) { + stream << curline << std::endl; + curline = std::string(indent_size, ' '); + } + curline += " ["; + } else if ((cur_mutex != nullptr) && (arg_mutex != nullptr)) { + if (cur_mutex != arg_mutex) { + curline += ']'; + if (this->m_usage_break_on_mutex || + curline.size() + 3 + arg_inline_usage.size() > + this->m_usage_max_line_width) { + stream << curline << std::endl; + curline = std::string(indent_size, ' '); + } + curline += " ["; + } else { + curline += '|'; + } + } + cur_mutex = arg_mutex; + if (curline.size() != indent_size && + curline.size() + 1 + arg_inline_usage.size() > + this->m_usage_max_line_width) { + stream << curline << std::endl; + curline = std::string(indent_size, ' '); + curline += " "; + } else if (cur_mutex == nullptr) { + curline += " "; + } + curline += arg_inline_usage; + } + if (cur_mutex != nullptr) { + curline += ']'; + } + return found_options; + }; + + const bool found_options = deal_with_options_of_group(0); + + if (found_options && multiline_usage && + !this->m_positional_arguments.empty()) { + stream << curline << std::endl; + curline = std::string(indent_size, ' '); + } + // Put positional arguments after the optionals + for (const auto &argument : this->m_positional_arguments) { + if (argument.m_is_hidden) { + continue; + } + const std::string pos_arg = !argument.m_metavar.empty() + ? argument.m_metavar + : argument.m_names.front(); + if (curline.size() + 1 + pos_arg.size() > this->m_usage_max_line_width) { + stream << curline << std::endl; + curline = std::string(indent_size, ' '); + } + curline += " "; + if (argument.m_num_args_range.get_min() == 0 && + !argument.m_num_args_range.is_right_bounded()) { + curline += "["; + curline += pos_arg; + curline += "]..."; + } else if (argument.m_num_args_range.get_min() == 1 && + !argument.m_num_args_range.is_right_bounded()) { + curline += pos_arg; + curline += "..."; + } else { + curline += pos_arg; + } + } + + if (multiline_usage) { + // Display options of other groups + for (std::size_t i = 0; i < m_group_names.size(); ++i) { + stream << curline << std::endl << std::endl; + stream << m_group_names[i] << ":" << std::endl; + curline = std::string(indent_size, ' '); + deal_with_options_of_group(i + 1); + } + } + + stream << curline; + + // Put subcommands after positional arguments + if (!m_subparser_map.empty()) { + stream << " {"; + std::size_t i{0}; + for (const auto &[command, subparser] : m_subparser_map) { + if (subparser->get().m_suppress) { + continue; + } + + if (i == 0) { + stream << command; + } else { + stream << "," << command; + } + ++i; + } + stream << "}"; + } + + return stream.str(); + } + + // Printing the one and only help message + // I've stuck with a simple message format, nothing fancy. + [[deprecated("Use cout << program; instead. See also help().")]] std::string + print_help() const { + auto out = help(); + std::cout << out.rdbuf(); + return out.str(); + } + + void add_subparser(ArgumentParser &parser) { + parser.m_parser_path = m_program_name + " " + parser.m_program_name; + auto it = m_subparsers.emplace(std::cend(m_subparsers), parser); + m_subparser_map.insert_or_assign(parser.m_program_name, it); + m_subparser_used.insert_or_assign(parser.m_program_name, false); + } + + void set_suppress(bool suppress) { m_suppress = suppress; } + +protected: + const MutuallyExclusiveGroup *get_belonging_mutex(const Argument *arg) const { + for (const auto &mutex : m_mutually_exclusive_groups) { + if (std::find(mutex.m_elements.begin(), mutex.m_elements.end(), arg) != + mutex.m_elements.end()) { + return &mutex; + } + } + return nullptr; + } + + bool is_valid_prefix_char(char c) const { + return m_prefix_chars.find(c) != std::string::npos; + } + + char get_any_valid_prefix_char() const { return m_prefix_chars[0]; } + + /* + * Pre-process this argument list. Anything starting with "--", that + * contains an =, where the prefix before the = has an entry in the + * options table, should be split. + */ + std::vector + preprocess_arguments(const std::vector &raw_arguments) const { + std::vector arguments{}; + for (const auto &arg : raw_arguments) { + + const auto argument_starts_with_prefix_chars = + [this](const std::string &a) -> bool { + if (!a.empty()) { + + const auto legal_prefix = [this](char c) -> bool { + return m_prefix_chars.find(c) != std::string::npos; + }; + + // Windows-style + // if '/' is a legal prefix char + // then allow single '/' followed by argument name, followed by an + // assign char, e.g., ':' e.g., 'test.exe /A:Foo' + const auto windows_style = legal_prefix('/'); + + if (windows_style) { + if (legal_prefix(a[0])) { + return true; + } + } else { + // Slash '/' is not a legal prefix char + // For all other characters, only support long arguments + // i.e., the argument must start with 2 prefix chars, e.g, + // '--foo' e,g, './test --foo=Bar -DARG=yes' + if (a.size() > 1) { + return (legal_prefix(a[0]) && legal_prefix(a[1])); + } + } + } + return false; + }; + + // Check that: + // - We don't have an argument named exactly this + // - The argument starts with a prefix char, e.g., "--" + // - The argument contains an assign char, e.g., "=" + auto assign_char_pos = arg.find_first_of(m_assign_chars); + + if (m_argument_map.find(arg) == m_argument_map.end() && + argument_starts_with_prefix_chars(arg) && + assign_char_pos != std::string::npos) { + // Get the name of the potential option, and check it exists + std::string opt_name = arg.substr(0, assign_char_pos); + if (m_argument_map.find(opt_name) != m_argument_map.end()) { + // This is the name of an option! Split it into two parts + arguments.push_back(std::move(opt_name)); + arguments.push_back(arg.substr(assign_char_pos + 1)); + continue; + } + } + // If we've fallen through to here, then it's a standard argument + arguments.push_back(arg); + } + return arguments; + } + + /* + * @throws std::runtime_error in case of any invalid argument + */ + void parse_args_internal(const std::vector &raw_arguments) { + auto arguments = preprocess_arguments(raw_arguments); + if (m_program_name.empty() && !arguments.empty()) { + m_program_name = arguments.front(); + } + auto end = std::end(arguments); + auto positional_argument_it = std::begin(m_positional_arguments); + for (auto it = std::next(std::begin(arguments)); it != end;) { + const auto ¤t_argument = *it; + if (Argument::is_positional(current_argument, m_prefix_chars)) { + if (positional_argument_it == std::end(m_positional_arguments)) { + + // Check sub-parsers + auto subparser_it = m_subparser_map.find(current_argument); + if (subparser_it != m_subparser_map.end()) { + + // build list of remaining args + const auto unprocessed_arguments = + std::vector(it, end); + + // invoke subparser + m_is_parsed = true; + m_subparser_used[current_argument] = true; + return subparser_it->second->get().parse_args( + unprocessed_arguments); + } + + if (m_positional_arguments.empty()) { + + // Ask the user if they argument they provided was a typo + // for some sub-parser, + // e.g., user provided `git totes` instead of `git notes` + if (!m_subparser_map.empty()) { + throw std::runtime_error( + "Failed to parse '" + current_argument + "', did you mean '" + + std::string{details::get_most_similar_string( + m_subparser_map, current_argument)} + + "'"); + } + + // Ask the user if they meant to use a specific optional argument + if (!m_optional_arguments.empty()) { + for (const auto &opt : m_optional_arguments) { + if (!opt.m_implicit_value.has_value()) { + // not a flag, requires a value + if (!opt.m_is_used) { + throw std::runtime_error( + "Zero positional arguments expected, did you mean " + + opt.get_usage_full()); + } + } + } + + throw std::runtime_error("Zero positional arguments expected"); + } else { + throw std::runtime_error("Zero positional arguments expected"); + } + } else { + throw std::runtime_error("Maximum number of positional arguments " + "exceeded, failed to parse '" + + current_argument + "'"); + } + } + auto argument = positional_argument_it++; + + // Deal with the situation of ... + if (argument->m_num_args_range.get_min() == 1 && + argument->m_num_args_range.get_max() == (std::numeric_limits::max)() && + positional_argument_it != std::end(m_positional_arguments) && + std::next(positional_argument_it) == std::end(m_positional_arguments) && + positional_argument_it->m_num_args_range.get_min() == 1 && + positional_argument_it->m_num_args_range.get_max() == 1 ) { + if (std::next(it) != end) { + positional_argument_it->consume(std::prev(end), end); + end = std::prev(end); + } else { + throw std::runtime_error("Missing " + positional_argument_it->m_names.front()); + } + } + + it = argument->consume(it, end); + continue; + } + + auto arg_map_it = m_argument_map.find(current_argument); + if (arg_map_it != m_argument_map.end()) { + auto argument = arg_map_it->second; + it = argument->consume(std::next(it), end, arg_map_it->first); + } else if (const auto &compound_arg = current_argument; + compound_arg.size() > 1 && + is_valid_prefix_char(compound_arg[0]) && + !is_valid_prefix_char(compound_arg[1])) { + ++it; + for (std::size_t j = 1; j < compound_arg.size(); j++) { + auto hypothetical_arg = std::string{'-', compound_arg[j]}; + auto arg_map_it2 = m_argument_map.find(hypothetical_arg); + if (arg_map_it2 != m_argument_map.end()) { + auto argument = arg_map_it2->second; + it = argument->consume(it, end, arg_map_it2->first); + } else { + throw std::runtime_error("Unknown argument: " + current_argument); + } + } + } else { + throw std::runtime_error("Unknown argument: " + current_argument); + } + } + m_is_parsed = true; + } + + /* + * Like parse_args_internal but collects unused args into a vector + */ + std::vector + parse_known_args_internal(const std::vector &raw_arguments) { + auto arguments = preprocess_arguments(raw_arguments); + + std::vector unknown_arguments{}; + + if (m_program_name.empty() && !arguments.empty()) { + m_program_name = arguments.front(); + } + auto end = std::end(arguments); + auto positional_argument_it = std::begin(m_positional_arguments); + for (auto it = std::next(std::begin(arguments)); it != end;) { + const auto ¤t_argument = *it; + if (Argument::is_positional(current_argument, m_prefix_chars)) { + if (positional_argument_it == std::end(m_positional_arguments)) { + + // Check sub-parsers + auto subparser_it = m_subparser_map.find(current_argument); + if (subparser_it != m_subparser_map.end()) { + + // build list of remaining args + const auto unprocessed_arguments = + std::vector(it, end); + + // invoke subparser + m_is_parsed = true; + m_subparser_used[current_argument] = true; + return subparser_it->second->get().parse_known_args_internal( + unprocessed_arguments); + } + + // save current argument as unknown and go to next argument + unknown_arguments.push_back(current_argument); + ++it; + } else { + // current argument is the value of a positional argument + // consume it + auto argument = positional_argument_it++; + it = argument->consume(it, end); + } + continue; + } + + auto arg_map_it = m_argument_map.find(current_argument); + if (arg_map_it != m_argument_map.end()) { + auto argument = arg_map_it->second; + it = argument->consume(std::next(it), end, arg_map_it->first); + } else if (const auto &compound_arg = current_argument; + compound_arg.size() > 1 && + is_valid_prefix_char(compound_arg[0]) && + !is_valid_prefix_char(compound_arg[1])) { + ++it; + for (std::size_t j = 1; j < compound_arg.size(); j++) { + auto hypothetical_arg = std::string{'-', compound_arg[j]}; + auto arg_map_it2 = m_argument_map.find(hypothetical_arg); + if (arg_map_it2 != m_argument_map.end()) { + auto argument = arg_map_it2->second; + it = argument->consume(it, end, arg_map_it2->first); + } else { + unknown_arguments.push_back(current_argument); + break; + } + } + } else { + // current argument is an optional-like argument that is unknown + // save it and move to next argument + unknown_arguments.push_back(current_argument); + ++it; + } + } + m_is_parsed = true; + return unknown_arguments; + } + + // Used by print_help. + std::size_t get_length_of_longest_argument() const { + if (m_argument_map.empty()) { + return 0; + } + std::size_t max_size = 0; + for ([[maybe_unused]] const auto &[unused, argument] : m_argument_map) { + max_size = + std::max(max_size, argument->get_arguments_length()); + } + for ([[maybe_unused]] const auto &[command, unused] : m_subparser_map) { + max_size = std::max(max_size, command.size()); + } + return max_size; + } + + using argument_it = std::list::iterator; + using mutex_group_it = std::vector::iterator; + using argument_parser_it = + std::list>::iterator; + + void index_argument(argument_it it) { + for (const auto &name : std::as_const(it->m_names)) { + m_argument_map.insert_or_assign(name, it); + } + } + + std::string m_program_name; + std::string m_version; + std::string m_description; + std::string m_epilog; + bool m_exit_on_default_arguments = true; + std::string m_prefix_chars{"-"}; + std::string m_assign_chars{"="}; + bool m_is_parsed = false; + std::list m_positional_arguments; + std::list m_optional_arguments; + std::map m_argument_map; + std::string m_parser_path; + std::list> m_subparsers; + std::map m_subparser_map; + std::map m_subparser_used; + std::vector m_mutually_exclusive_groups; + bool m_suppress = false; + std::size_t m_usage_max_line_width = (std::numeric_limits::max)(); + bool m_usage_break_on_mutex = false; + int m_usage_newline_counter = 0; + std::vector m_group_names; +}; + +} // namespace argparse diff --git a/src/project3/include/classes.hpp b/src/project3/include/classes.hpp new file mode 100644 index 0000000..064093d --- /dev/null +++ b/src/project3/include/classes.hpp @@ -0,0 +1,85 @@ +#include +#include + +#ifndef CLASSES_HPP +#define CLASSES_HPP + +class PenningTrap; + +arma::vec standard_external_E_field(const arma::vec& r, double t, const PenningTrap& trap); +arma::vec oscillating_external_E_field(const arma::vec& r, double t, const PenningTrap& trap); +arma::vec standard_external_B_field(const arma::vec& r, double t, const PenningTrap& trap); +arma::vec limited_range_external_B_field(const arma::vec& r, double t, const PenningTrap& trap); + +class Particle { + public: + arma::vec position; + arma::vec velocity; + double mass; + double charge; + Particle(arma::vec pos, arma::vec vel, double m, double q); + Particle(double m, double q); + std::string info() const; + bool inside_trap(double d) const; + +}; + +class PenningTrap { + private: + std::vector particles; + bool coulomb_interaction = true; + bool modified_fields = false; + arma::vec (*_external_E_field)(const arma::vec& r, double t, const PenningTrap& trap) = standard_external_E_field; + arma::vec (*_external_B_field)(const arma::vec& r, double t, const PenningTrap& trap) = standard_external_B_field; + public: + double B0; + double V0; + double d; + double f = 0.0; // Amplitude of oscillating potential + double omega_V = 0.0; // Angular frequency of oscillating potential + PenningTrap(double B0_in, double V0_in, double d_in); + void add_particle(Particle& p); + void add_n_identical_particles(int n, double m, double q); + arma::vec external_E_field(const arma::vec& r, double t) const { + return _external_E_field(r, t, *this); + } + arma::vec external_B_field(const arma::vec& r, double t) const { + return _external_B_field(r, t, *this); + } + void set_external_E_field(arma::vec (*E_field_func)(const arma::vec&, double, const PenningTrap&)); + void set_external_B_field(arma::vec (*B_field_func)(const arma::vec&, double, const PenningTrap&)); + arma::vec force_particle(int i, int j) const; + arma::vec total_force_external(int i, double t) const; + arma::vec total_force_particles(int i) const; + arma::vec total_force(int i, double t) const; + double potential_energy(int i) const; + double total_potential_energy() const; + double kinetic_energy(int i) const; + double total_kinetic_energy() const; + double total_energy() const; + Particle& get_particle(int i); + const Particle& get_particle(int i) const; + int num_particles() const; + void enable_interactions() { + coulomb_interaction = true; + } + void disable_interactions() { + coulomb_interaction = false; + } + bool get_coulomb_interaction() const { + return coulomb_interaction; + } + bool get_modified_fields() const { + return modified_fields; + } + std::string info() const; + int num_inside_trap() const; +}; + + +std::ostream& operator<<(std::ostream& os, const Particle& p); + +std::ostream& operator<<(std::ostream& os, const PenningTrap& trap); + + +#endif \ No newline at end of file diff --git a/src/project3/include/constants.hpp b/src/project3/include/constants.hpp new file mode 100644 index 0000000..55fd7fb --- /dev/null +++ b/src/project3/include/constants.hpp @@ -0,0 +1,13 @@ +#ifndef CONSTANTS_HPP +#define CONSTANTS_HPP + +namespace constants { + const double k_e = 1.38935333e5; // Unit: u * um^3 / (us^2 * e^2), Coulomb's constant + const double tesla = 9.64852558e1; // Conversion factor from tesla to u/(us * e) + const double volt = 9.64852558e7; // Conversion factor from volt to u * um^2 / (us^2 * e) + const double meter = 1e6; // Conversion factor from meter to micrometer + const double second = 1e6; // Conversion factor from second to microsecond + const double amu = 1.0; // Atomic mass unit in terms of itself (u) + const double elementary_charge = 1.0; // Elementary charge in terms of itself (e) +} +#endif \ No newline at end of file diff --git a/src/project3/include/solvers.hpp b/src/project3/include/solvers.hpp new file mode 100644 index 0000000..32c11cc --- /dev/null +++ b/src/project3/include/solvers.hpp @@ -0,0 +1,92 @@ +#include "classes.hpp" +#include +#include + +#ifndef SOLVERS_HPP +#define SOLVERS_HPP + + +arma::cube to_cube(const std::vector& vec); + +void save_to_file(const std::string& filename, const arma::cube& data); +void save_to_file(const std::string& filename, const arma::vec& data); + + +class Solver { + protected: + PenningTrap& trap; + double dt; + int steps; + std::vector positions; + std::vector velocities; + std::vector inside_trap_count; + std::vector potential_energies; + std::vector kinetic_energies; + std::vector total_energies; + virtual void step() = 0; + void record_position(); + void record_velocity(); + void record_energies(); + void record_trap_count(); + public: + Solver(PenningTrap& trap_in, double dt_in); + void simulate(int num_steps); + void simulate_time(double total_time) { + int num_steps = total_steps(total_time); + simulate(num_steps); + } + std::vector get_positions(); + std::vector get_velocities(); + std::vector get_potential_energies(); + std::vector get_kinetic_energies(); + std::vector get_total_energies(); + std::vector get_inside_trap_count(); + int get_current_step(); + int total_steps(double total_time); + void save(std::string base_filename); + void save_counts(std::string base_filename); + void record_current_step() { + record_position(); + record_velocity(); + if (!trap.get_modified_fields()) { + record_energies(); + } + record_trap_count(); + } + bool record_positions_and_velocities = true; + bool record_inside_trap_count = true; +}; + +class EulerSolver : public Solver { + public: + EulerSolver(PenningTrap& trap_in, double dt_in); + void step() override; +}; +class RK4Solver : public Solver { + public: + RK4Solver(PenningTrap& trap_in, double dt_in); + void step() override; +}; +class VelocityVerletSolver : public Solver { + private: + std::vector accelerations; + public: + VelocityVerletSolver(PenningTrap& trap_in, double dt_in); + void step() override; +}; + +class AnalyticalSolver : public Solver { + private: + std::vector w_0; + std::vector w_z; + std::vector w_plus; + std::vector w_minus; + std::vector A_plus; + std::vector A_minus; + std::vector z_0; + public: + AnalyticalSolver(PenningTrap& trap_in, double dt_in); + void step() override; +}; + +#endif \ No newline at end of file diff --git a/src/project3/many_particles.cpp b/src/project3/many_particles.cpp new file mode 100644 index 0000000..19efa67 --- /dev/null +++ b/src/project3/many_particles.cpp @@ -0,0 +1,173 @@ +#include +#include "classes.hpp" +#include "solvers.hpp" +#include "constants.hpp" +#include "argparse/argparse.hpp" + +#include +#include + +using namespace std; + +void print_first_particles(PenningTrap& trap, int num_particles) { + for (int i = 0; i < num_particles && i < trap.num_particles(); ++i) { + cout << "Particle " << i + 1 << ": " << endl << trap.get_particle(i) << endl; + } +} + +void save_timing(clock_t start, clock_t end, const string& filename) { + double cpu_time_used = ((double)(end - start)) / CLOCKS_PER_SEC; + ofstream file(filename); + if (file.is_open()) { + file << "CPU time used: " << cpu_time_used << " seconds" << endl; + file.close(); + } else { + cerr << "Unable to open file for writing timing information." << endl; + } +} + +int main(int argc, char* argv[]) { + argparse::ArgumentParser program("many_particles"); + program.add_argument("-v", "--verbose") + .help("Enable verbose output of first particles") + .flag(); + program.add_argument("-N", "--steps") + .help("Number of time steps") + .default_value(4000) + .scan<'i', int>(); + program.add_argument("-n", "--num_particles") + .help("Number of particles") + .default_value(100) + .scan<'i', int>(); + program.add_argument("-i", "--disable-interactions") + .help("Disable Coulomb interactions") + .default_value(true) + .implicit_value(false); + program.add_argument("-O", "--oscillating-potential") + .help("Enable oscillating potential") + .flag(); + program.add_argument("-f", "--amplitude") + .help("Amplitude of oscillating potential") + .default_value(0.0) + .scan<'f', double>(); + program.add_argument("-w", "--frequency") + .help("Angular frequency of oscillating potential") + .default_value(0.0) + .scan<'f', double>(); + program.add_argument("-r", "--reduced-output") + .help("Reduce output for oscillating potential to only final count") + .flag(); + program.add_argument("-E", "--euler") + .help("Use Euler solver instead of RK4") + .flag(); + program.add_argument("-V", "--velocity-verlet") + .help("Use Velocity Verlet solver instead of RK4") + .flag(); + program.add_argument("-t", "--timing") + .help("Enable timing of the simulation") + .flag(); + + try { + program.parse_args(argc, argv); + } catch (const std::runtime_error& err) { + std::cerr << err.what() << std::endl; + std::cerr << program.help().str() << std::endl; + exit(1); + } + + bool verbose = program.get("--verbose"); + + int N = program.get("--steps"); + int n_particles = program.get("--num_particles"); + bool interactions = program.get("--disable-interactions"); + bool oscillating_potential = program.get("--oscillating-potential"); + bool use_euler = program.get("--euler"); + bool use_velocity_verlet = program.get("--velocity-verlet"); + bool timing = program.get("--timing"); + + double B_0 = 1.0; // Tesla + double V_0 = 0.025; // Volt + double d = 500e-6; // meter + + PenningTrap trap(B_0, V_0, d); + if (interactions) { + trap.enable_interactions(); + } + else { + trap.disable_interactions(); + } + + + bool reduced_output = program.get("--reduced-output"); // Only relevant if oscillating potential is enabled + if (oscillating_potential) { + trap.set_external_E_field(oscillating_external_E_field); + trap.set_external_B_field(limited_range_external_B_field); + trap.f = program.get("--amplitude"); + trap.omega_V = program.get("--frequency"); + } + + trap.add_n_identical_particles(n_particles, 40.078 * constants::amu, 1 * constants::elementary_charge); + if (verbose) { + cout << trap.info() << endl; + print_first_particles(trap, 10); + } + + + double dt = 50e-6 / N; // seconds + + + unique_ptr solver; + string solver_name; + + if (use_euler) { + solver = make_unique(trap, dt); + solver_name = "Euler"; + } else if (use_velocity_verlet) { + solver = make_unique(trap, dt); + solver_name = "VelocityVerlet"; + } else { + solver = make_unique(trap, dt); + solver_name = "RK4"; + } + clock_t start, end; + start = clock(); + if (oscillating_potential) { + solver->record_positions_and_velocities = false; // Disable recording for performance + if (reduced_output) { + solver->record_inside_trap_count = false; // Disable recording for performance + } + solver->simulate(N); + end = clock(); + string folder = "results/" + solver_name + "/many_particles_" + to_string(trap.num_particles()) + "_particles_" + to_string(int(N)) + "_steps_" + (interactions ? "enabled" : "disabled") + "_interactions_oscillating_potential"; + string filename = "f" + to_string(trap.f) + "_w" + to_string(trap.omega_V); + filesystem::create_directories(folder); + + if (reduced_output) { + solver->record_current_step(); // Record only final step + } + solver->save_counts(folder + "/" + filename); + cout << "Final number of particles inside trap: " << trap.num_inside_trap() << " out of " << trap.num_particles() << endl; + if (timing) { + save_timing(start, end, folder + "/" + filename + "_timing.txt"); + } + return 0; + + } + solver->simulate(N); + end = clock(); + string folder = "results/" + solver_name; + filesystem::create_directories(folder); + string filename = "/many_particles_" + to_string(trap.num_particles()) + "_particles_" + to_string(int(N)) + "_steps_" + (interactions ? "enabled" : "disabled") + "_interactions"; + solver->save(folder + filename); + if (verbose) { + vector positions = solver->get_positions(); + cout << "Final positions of particles: " << endl; + print_first_particles(trap, 10); + } + if (timing) { + save_timing(start, end, folder + filename + "_timing.txt"); + } + + return 0; +} + diff --git a/src/project3/progress.sh b/src/project3/progress.sh new file mode 100755 index 0000000..e9155df --- /dev/null +++ b/src/project3/progress.sh @@ -0,0 +1,36 @@ +#!/usr/bin/env bash + +interval=5 # seconds between updates + +while true; do + data=$(pueue status --json) + + queued=$(jq '.tasks | map(select(.status=="Queued")) | length' <<< "$data") + running=$(jq '.tasks | map(select(.status=="Running")) | length' <<< "$data") + total=$(jq '.tasks | length' <<< "$data") + + # average duration of finished tasks (seconds) + avg=$(jq '.tasks + | map(select(.start != null and .end != null) + | ((.end[0:19] + "Z") | fromdate) - ((.start[0:19] + "Z") | fromdate)) + | if length > 0 then add / length else 0 end' <<< "$data") + + clear + echo "Queued: $queued" + echo "Running: $running" + echo "Total: $total" + + # Only calculate ETA if avg > 0 + if (( $(echo "$avg > 0" | bc -l) )); then + remaining=$((queued + running)) + eta_seconds=$(echo "$remaining * $avg" | bc -l) + finish_time=$(date -d "@$(($(date +%s) + ${eta_seconds%.*}))" +"%Y-%m-%d %H:%M:%S") + printf "Average duration: %.2fs\n" "$avg" + printf "ETA: ~%.0fs (finishing around %s)\n" "$eta_seconds" "$finish_time" + else + echo "ETA: Not enough finished jobs in history yet." + fi + + sleep $interval +done + diff --git a/src/project3/python/plot_energies.py b/src/project3/python/plot_energies.py new file mode 100644 index 0000000..0b7b959 --- /dev/null +++ b/src/project3/python/plot_energies.py @@ -0,0 +1,35 @@ +import matplotlib.pyplot as plt +import numpy as np +import os +import typer + +app = typer.Typer() + +def load_data(total_energy_file): + total_energy = np.loadtxt(total_energy_file) + kin_file = total_energy_file.replace("total_energies", "kinetic_energies") + pot_file = total_energy_file.replace("total_energies", "potential_energies") + kinetic_energy = np.loadtxt(kin_file) + potential_energy = np.loadtxt(pot_file) + step_count = len(total_energy) + return total_energy, kinetic_energy, potential_energy, step_count + +@app.command() +def plot_kinetic_energy(energy_files: list[str]=[], labels: list[str]=[], mass: float=1.0): + """Plot kinetic energy from velocity files.""" + for i, energy_file in enumerate(energy_files): + total_energy, kinetic_energy, potential_energy, step_count = load_data(energy_file) + time = np.linspace(0, 50e-6, step_count) + color = f"C{i % 10}" + plt.plot(time, total_energy, label=labels[i], color=color, linestyle="-") + plt.plot(time, potential_energy, color=color, linestyle="--") + plt.plot(time, kinetic_energy, color=color, linestyle=":") + plt.xlabel('Time (s)') + plt.ylabel('Energy (J)') + plt.title('Energy Evolution') + plt.legend() + plt.show() + + +if __name__ == "__main__": + app() \ No newline at end of file diff --git a/src/project3/python/plot_resonances.py b/src/project3/python/plot_resonances.py new file mode 100644 index 0000000..359acbb --- /dev/null +++ b/src/project3/python/plot_resonances.py @@ -0,0 +1,49 @@ +import matplotlib.pyplot as plt +import numpy as np +import os +import argparse + +argparser = argparse.ArgumentParser() +argparser.add_argument("directory", type=str, help="Directory containing the data files") +args = argparser.parse_args() +directory = args.directory + +def get_filenames(directory): + filenames = [os.path.join(directory, f) for f in os.listdir(directory) if f.endswith(".red")] + return filenames + +def get_pars(filename): + amplitude = float(filename.split("_")[0][1:]) + omega = float(filename.split("_")[1][1:]) + return amplitude, omega + +def load_data(filename): + with open(filename, "r") as file: + x = int(file.readline().split()[1]) + return x + +def plot_resonances(amplitudes, omegas, counts): + plt.figure(figsize=(8, 6)) + unique_amps = sorted(set(amplitudes)) + for amp in unique_amps: + mask = np.array(amplitudes) == amp + omegas_subset = np.array(omegas)[mask] + counts_subset = np.array(counts)[mask] + omegas_sorted, counts_sorted = zip(*sorted(zip(omegas_subset, counts_subset))) + plt.plot(omegas_sorted, counts_sorted, label=f'Amplitude {amp}') + plt.xlabel('Omega') + plt.ylabel('Count inside trap') + plt.legend() + plt.show() + +if __name__ == "__main__": + filenames = get_filenames(directory) + amplitudes, omegas, counts = [], [], [] + for filename in filenames: + amplitude, omega = get_pars(os.path.basename(filename)) + count = load_data(filename) + amplitudes.append(amplitude) + omegas.append(omega) + counts.append(count) + plot_resonances(amplitudes, omegas, counts) + diff --git a/src/project3/python/plot_trajectories.py b/src/project3/python/plot_trajectories.py new file mode 100644 index 0000000..620924d --- /dev/null +++ b/src/project3/python/plot_trajectories.py @@ -0,0 +1,93 @@ +import matplotlib.pyplot as plt +import numpy as np +import os +import typer +from typing import Literal + +app = typer.Typer() + +def load_trajectory_data(filename): + data = np.loadtxt(filename) + step_count = int(float(filename.split("/")[-1].split("_")[5])) + particle_count = int(float(filename.split("/")[-1].split("_")[2])) + single_particle = particle_count == 1 + if single_particle: + data = data.reshape((data.shape[0], 1)) + data = data.reshape(step_count, 3, particle_count).transpose(1, 2, 0) + return data, step_count, particle_count + +@app.command() +def plot_z_evolution(trajectory_file: str): + """Plot the z-coordinate evolution from a trajectory file.""" + data, step_count, particle_count = load_trajectory_data(trajectory_file) + time = np.linspace(0, 50e-6, step_count) + for i in range(particle_count): + plt.plot(time, data[2, i, :] * 1e-6, label=f'Particle {i+1}') + plt.xlabel('Time (s)') + plt.ylabel('Z Position (m)') + plt.title('Z Coordinate Evolution') + plt.legend() + plt.show() + +@app.command() +def plot_xy_trajectory(trajectory_file: str, plot_start_end: bool = True, plot_time: bool = False): + """Plot the XY trajectory from a trajectory file.""" + data, step_count, particle_count = load_trajectory_data(trajectory_file) + for i in range(particle_count): + if plot_start_end: + plt.scatter(data[0, i, 0] * 1e-6, data[1, i, 0] * 1e-6, marker="o", color="black", label="Start" if i == 0 else "") + plt.scatter(data[0, i, -1] * 1e-6, data[1, i, -1] * 1e-6, marker="x", color="red", label="End" if i == 0 else "") + if plot_time: + marker = [".", ",", "o", "v", "^", "<", ">", "1", "2", "3", "4", "8"][i % 12] + cm = plt.get_cmap('viridis') + c = cm(np.linspace(0, 1, step_count)) + plt.scatter(data[0, i, :] * 1e-6, data[1, i, :] * 1e-6, label=f'Particle {i+1}', marker=marker, c=c) + else: + ls = "-" + c = f"C{i % 10}" + plt.plot(data[0, i, :] * 1e-6, data[1, i, :] * 1e-6, label=f'Particle {i+1}', linestyle=ls, color=c) + plt.xlabel('X Position (m)') + plt.ylabel('Y Position (m)') + plt.title('XY Trajectory') + plt.axis('equal') + if plot_time: + #sm = plt.cm.ScalarMappable(cmap=cm, norm=plt.Normalize(vmin=0, vmax=50e-6)) + #sm.set_array([]) + #cbar = plt.colorbar(sm) + #cbar.set_label('Time') + pass + plt.legend() + plt.show() + +@app.command() +def phase_space_plot(trajectory_file: str, velocity_file: str, axis: Literal['x', 'y', 'z'] = 'z'): + """Plot phase space (Z vs Vz) from trajectory and velocity files.""" + axis_index = {'x': 0, 'y': 1, 'z': 2}[axis] + traj_data, step_count, particle_count = load_trajectory_data(trajectory_file) + vel_data, _, _ = load_trajectory_data(velocity_file) + for i in range(particle_count): + plt.plot(traj_data[axis_index, i, :] * 1e-6, vel_data[axis_index, i, :] * 1e-6, label=f'Particle {i+1}') + plt.xlabel(f'{axis.upper()} Position (m)') + plt.ylabel(f'{axis.upper()} Velocity (m/s)') + plt.title(f'Phase Space Plot ({axis.upper()} vs V{axis.upper()})') + plt.legend() + plt.show() + +@app.command() +def plot_error(trajectory_files: list[str], analytical_file: str): + """Plot error between numerical and analytical trajectories.""" + analytical_data, step_count, particle_count = load_trajectory_data(analytical_file) + time = np.linspace(0, 50e-6, step_count) + for traj_file in trajectory_files: + numerical_data, _, _ = load_trajectory_data(traj_file) + error = numerical_data - analytical_data + relative_error = np.mean(np.linalg.norm(error, axis=0) / np.linalg.norm(analytical_data, axis=0), axis=0) + plt.plot(time, relative_error, label=os.path.basename(traj_file)) + plt.xlabel('Time (s)') + plt.ylabel('Relative Error') + plt.title('Relative Error between Numerical and Analytical Trajectories') + plt.legend() + plt.show() + +if __name__ == "__main__": + app() \ No newline at end of file diff --git a/src/project3/src/classes.cpp b/src/project3/src/classes.cpp new file mode 100644 index 0000000..da678f4 --- /dev/null +++ b/src/project3/src/classes.cpp @@ -0,0 +1,223 @@ +#include "classes.hpp" +#include "constants.hpp" +#include + +// Initialize particle with position and velocity +// Mass in atomic mass units (u), charge in elementary charges (e) +// Positions in meters, velocities in meters/second +Particle::Particle(arma::vec pos, arma::vec vel, double m, double q) + : position(pos * constants::meter), velocity(vel * constants::meter / constants::second), mass(m * constants::amu), charge(q * constants::elementary_charge) {} + +// Initialize particle at origin with zero velocity +Particle::Particle(double m, double q) + : position(arma::vec({0,0,0})), velocity(arma::vec({0,0,0})), mass(m * constants::amu), charge(q * constants::elementary_charge) {} + +std::string Particle::info() const { + std::string info_str = "Particle (m = " + std::to_string(mass) + " u, q = " + std::to_string(charge) + " e)\n"; + info_str += "Position (um): [" + std::to_string(position(0)) + ", " + std::to_string(position(1)) + ", " + std::to_string(position(2)) + "]\n"; + info_str += "Velocity (um/us): [" + std::to_string(velocity(0)) + ", " + std::to_string(velocity(1)) + ", " + std::to_string(velocity(2)) + "]\n"; + return info_str; +} + +bool Particle::inside_trap(double d) const { + return arma::norm(position) <= d; +} + +// Empty Penning Trap with specified field parameters +// B0 in tesla, V0 in volt, d in meters +PenningTrap::PenningTrap(double B0_in, double V0_in, double d_in) + : B0(B0_in * constants::tesla), V0(V0_in * constants::volt), d(d_in * constants::meter) {} + +// Add particle to trap +void PenningTrap::add_particle(Particle& p) { + particles.push_back(p); +} + +// Add n identical particles with mass m and charge q +void PenningTrap::add_n_identical_particles(int n, double m, double q) { + for (int i = 0; i < n; ++i) { + arma::vec pos = arma::vec(3).randn() * 0.1 * d; // Random position within 10% of d (unit: um) + arma::vec vel = arma::vec(3).randn() * 0.1 * d; // Random velocity within 10% of d/s (unit: um/us) + Particle p(pos / constants::meter, vel, m, q); + add_particle(p); + } +} + +void PenningTrap::set_external_E_field(arma::vec (*E_field_func)(const arma::vec&, double, const PenningTrap&)) { + _external_E_field = E_field_func; + modified_fields = true; +} + +void PenningTrap::set_external_B_field(arma::vec (*B_field_func)(const arma::vec&, double, const PenningTrap&)) { + _external_B_field = B_field_func; + modified_fields = true; +} + + +// Force on particle i from particle j due to Coulomb interaction +// Neglecting magnetic fields from moving charges +arma::vec PenningTrap::force_particle(int i, int j) const { + if (!coulomb_interaction) { + return arma::vec({0, 0, 0}); + } + arma::vec r_i = particles[i].position; + arma::vec r_j = particles[j].position; + arma::vec r_ij = r_i - r_j; + double distance = arma::norm(r_ij); + if (distance == 0) { + return arma::vec({0, 0, 0}); // Avoid division by zero + } + double force_magnitude_per_dist = (constants::k_e * particles[i].charge * particles[j].charge) / (distance * distance * distance); + return force_magnitude_per_dist * r_ij; +} + +// Total force on particle i from external fields +// Sum of electric and magnetic forces +arma::vec PenningTrap::total_force_external(int i, double t) const { + arma::vec E = external_E_field(get_particle(i).position, t); + arma::vec B = external_B_field(get_particle(i).position, t); + arma::vec v = get_particle(i).velocity; + arma::vec F_electric = get_particle(i).charge * E; + arma::vec F_magnetic = get_particle(i).charge * arma::cross(v, B); + return F_electric + F_magnetic; +} + +// Total force on particle i from all other particles +arma::vec PenningTrap::total_force_particles(int i) const { + arma::vec total_force = arma::vec({0, 0, 0}); + if (!coulomb_interaction) { // Add the check here as well to avoid unnecessary calls to force_particle + return total_force; + } + for (size_t j = 0; j < particles.size(); ++j) { + if (j != i) { + total_force += force_particle(i, j); + } + } + return total_force; +} + +// Total force on particle i from all other particles and external fields +arma::vec PenningTrap::total_force(int i, double t) const { + return total_force_external(i, t) + total_force_particles(i); +} + +// Potential energy of particle i in the trap +double PenningTrap::potential_energy(int i) const { + if (modified_fields) { + throw std::runtime_error("Potential energy calculation not valid with modified external fields."); + } + double q = particles[i].charge; + arma::vec r = particles[i].position; + double V = (V0 / (d * d)) * (r(2) * r(2) - 0.5 * (r(0) * r(0) + r(1) * r(1))); + if (coulomb_interaction) { + // Add contributions from other particles + for (size_t j = 0; j < particles.size(); ++j) { + if (j != i) { + arma::vec r_j = particles[j].position; + double distance = arma::norm(r - r_j); + if (distance != 0) { + V += (constants::k_e * particles[j].charge) / distance; + } + } + } + } + return q * V; +} + +// Total potential energy of all particles in the trap +double PenningTrap::total_potential_energy() const { + if (modified_fields) { + throw std::runtime_error("Potential energy calculation not valid with modified external fields."); + } + double total_U = 0.0; + for (size_t i = 0; i < particles.size(); ++i) { + total_U += potential_energy(i); + } + return total_U; +} + +// Kinetic energy of particle i +double PenningTrap::kinetic_energy(int i) const { + double m = particles[i].mass; + arma::vec v = particles[i].velocity; + return 0.5 * m * arma::dot(v, v); +} +// Total kinetic energy of all particles in the trap +double PenningTrap::total_kinetic_energy() const { + double total_K = 0.0; + for (size_t i = 0; i < particles.size(); ++i) { + total_K += kinetic_energy(i); + } + return total_K; +} + +// Total energy of all particles in the trap +double PenningTrap::total_energy() const { + return total_kinetic_energy() + total_potential_energy(); +} + + +// Getter function for particles +Particle& PenningTrap::get_particle(int i) { + return particles[i]; +} + +const Particle& PenningTrap::get_particle(int i) const { + return particles[i]; +} + +// Number of particles in trap +int PenningTrap::num_particles() const{ + return particles.size(); +} + + +std::string PenningTrap::info() const { + std::string info_str = "Penning Trap (B0 = " + std::to_string(B0) + " u/(us*e), V0 = " + std::to_string(V0) + " u*um^2/(us^2*e), d = " + std::to_string(d) + " um)\n"; + info_str += "Number of particles: " + std::to_string(num_particles()) + "\n"; + info_str += "Coulomb interactions: " + std::string(coulomb_interaction ? "enabled" : "disabled") + "\n"; + return info_str; +} + +int PenningTrap::num_inside_trap() const { + int count = 0; + for (int i = 0; i < num_particles(); ++i) { + if (particles[i].inside_trap(d)) { + count++; + } + } + return count; +} + +// Electric field at position r and time t +// Electric field from V = V0 * (z^2 - (x^2 + y^2)/2) / d^2 +arma::vec standard_external_E_field(const arma::vec& r, double t, const PenningTrap& trap) { + return (trap.V0 / (trap.d * trap.d)) * arma::vec({r(0), r(1), -2 * r(2)}); +} + +// Magnetic field at position r and time t +// Constant magnetic field in z-direction +arma::vec standard_external_B_field(const arma::vec& r, double t, const PenningTrap& trap) { + return arma::vec({0, 0, trap.B0}); +} + +arma::vec oscillating_external_E_field(const arma::vec& r, double t, const PenningTrap& trap) { + if (arma::norm(r) > trap.d) { + return arma::vec({0, 0, 0}); // No electric field outside the trap + } + return standard_external_E_field(r, t, trap) * (1 + trap.f * cos(trap.omega_V * t)); +} + +arma::vec limited_range_external_B_field(const arma::vec& r, double t, const PenningTrap& trap) { + if (arma::norm(r) > trap.d) { + return arma::vec({0, 0, 0}); // No magnetic field outside the trap + } + return standard_external_B_field(r, t, trap); +} + +std::ostream& operator<<(std::ostream& os, const Particle& p) { + return os << p.info(); +} +std::ostream& operator<<(std::ostream& os, const PenningTrap& trap) { + return os << trap.info(); +} \ No newline at end of file diff --git a/src/project3/src/solvers.cpp b/src/project3/src/solvers.cpp new file mode 100644 index 0000000..5a88cdc --- /dev/null +++ b/src/project3/src/solvers.cpp @@ -0,0 +1,340 @@ +#include "classes.hpp" +#include "solvers.hpp" +#include "constants.hpp" + + +// Initialize base solver class +// dt in seconds +Solver::Solver(PenningTrap& trap_in, double dt_in) + : trap(trap_in), dt(dt_in * constants::second), steps(0) {} + +// Simulate the system for a given number of steps +void Solver::simulate(int num_steps) { + for (int i = 0; i < num_steps; i++) { + step(); + if (record_positions_and_velocities) { // Allow disabling of recording for performance + // Store positions of all particles after each step + record_position(); + // Store velocities of all particles after each step + record_velocity(); + if (!trap.get_modified_fields()) { + record_energies(); + } + } + if (record_inside_trap_count) { + record_trap_count(); + } + steps++; + } +} + +void Solver::record_position() { + arma::mat current_positions(3, trap.num_particles()); + for (int j = 0; j < trap.num_particles(); j++) { + current_positions.col(j) = trap.get_particle(j).position; + } + positions.push_back(current_positions); +} +void Solver::record_velocity() { + arma::mat current_velocities(3, trap.num_particles()); + for (int j = 0; j < trap.num_particles(); j++) { + current_velocities.col(j) = trap.get_particle(j).velocity; + } + velocities.push_back(current_velocities); +} + +void Solver::record_energies() { + double U = trap.total_potential_energy(); + double K = trap.total_kinetic_energy(); + double E = trap.total_energy(); + potential_energies.push_back(U); + kinetic_energies.push_back(K); + total_energies.push_back(E); +} + +void Solver::record_trap_count() { + inside_trap_count.push_back(trap.num_inside_trap()); +} + +// Get recorded positions of particles +std::vector Solver::get_positions() { + return positions; +} + +// Get recorded velocities of particles +std::vector Solver::get_velocities() { + return velocities; +} +std::vector Solver::get_potential_energies() { + return potential_energies; +} +std::vector Solver::get_kinetic_energies() { + return kinetic_energies; +} +std::vector Solver::get_total_energies() { + return total_energies; +} + +std::vector Solver::get_inside_trap_count() { + return inside_trap_count; +} + +// Get the current step count +int Solver::get_current_step() { + return steps; +} + +// Calculate total steps for a given total time +// total_time in seconds +int Solver::total_steps(double total_time) { + return static_cast(total_time * constants::second / dt); +} + +void Solver::save(std::string base_filename) { + arma::cube pos_cube = to_cube(positions); + save_to_file(base_filename + "_positions.dat", pos_cube); + arma::cube vel_cube = to_cube(velocities); + save_to_file(base_filename + "_velocities.dat", vel_cube); + if (!trap.get_modified_fields()) { + arma::vec U = arma::vec(potential_energies); + save_to_file(base_filename + "_potential_energies.dat", U); + arma::vec K = arma::vec(kinetic_energies); + save_to_file(base_filename + "_kinetic_energies.dat", K); + arma::vec E = arma::vec(total_energies); + save_to_file(base_filename + "_total_energies.dat", E); + } +} + +void Solver::save_counts(std::string base_filename) { + std::string extension = ".dat"; + if (inside_trap_count.size() == 1) { + // Reduced to last element if only one entry + extension = ".red"; + } + std::ofstream file(base_filename + "_inside_trap_count" + extension); + for (size_t i = 0; i < inside_trap_count.size(); i++) { + file << i << " " << inside_trap_count[i] << "\n"; + } + file.close(); +} + +// Initialize Euler solver +EulerSolver::EulerSolver(PenningTrap& trap_in, double dt_in) + : Solver(trap_in, dt_in) {} + +// Perform a single Euler integration step +void EulerSolver::step() { + int n = trap.num_particles(); + for (int i = 0; i < n; i++) { + Particle& p = trap.get_particle(i); + arma::vec F = trap.total_force(i, steps * dt); + arma::vec a = F / p.mass; + p.position += p.velocity * dt; + p.velocity += a * dt; + } +} + +// Initialize Runge-Kutta 4 solver +RK4Solver::RK4Solver(PenningTrap& trap_in, double dt_in) + : Solver(trap_in, dt_in) {} + + +// Perform a single RK4 integration step +void RK4Solver::step() { + int n = trap.num_particles(); + std::vector k1_v(n), k1_r(n); + std::vector k2_v(n), k2_r(n); + std::vector k3_v(n), k3_r(n); + std::vector k4_v(n), k4_r(n); + + std::vector current_positions(n); + std::vector current_velocities(n); + + for (int i = 0; i < n; i++) { + current_positions[i] = trap.get_particle(i).position; + current_velocities[i] = trap.get_particle(i).velocity; + } + double t = steps * dt; + // Calculate k1 + for (int i = 0; i < n; i++) { + Particle& p = trap.get_particle(i); + arma::vec F = trap.total_force(i, t); + arma::vec a = F / p.mass; + k1_v[i] = a * dt; + k1_r[i] = p.velocity * dt; + p.position += k1_r[i] / 2; + p.velocity += k1_v[i] / 2; + } + // Calculate k2 + for (int i = 0; i < n; i++) { + Particle& p = trap.get_particle(i); + arma::vec F = trap.total_force(i, t + dt / 2); + arma::vec a = F / p.mass; + k2_v[i] = a * dt; + k2_r[i] = p.velocity * dt; + p.position = current_positions[i] + k2_r[i] / 2; + p.velocity = current_velocities[i] + k2_v[i] / 2; + } + // Calculate k3 + for (int i = 0; i < n; i++) { + Particle& p = trap.get_particle(i); + arma::vec F = trap.total_force(i, t + dt / 2); + arma::vec a = F / p.mass; + k3_v[i] = a * dt; + k3_r[i] = p.velocity * dt; + p.position = current_positions[i] + k3_r[i]; + p.velocity = current_velocities[i] + k3_v[i]; + } + // Calculate k4 + for (int i = 0; i < n; i++) { + Particle& p = trap.get_particle(i); + arma::vec F = trap.total_force(i, t + dt); + arma::vec a = F / p.mass; + k4_v[i] = a * dt; + k4_r[i] = p.velocity * dt; + p.position = current_positions[i]; + p.velocity = current_velocities[i]; + } + // Update positions and velocities + for (int i = 0; i < n; i++) { + Particle& p = trap.get_particle(i); + p.position = current_positions[i] + (k1_r[i] + 2.0 * k2_r[i] + 2.0 * k3_r[i] + k4_r[i]) / 6.0; + p.velocity = current_velocities[i] + (k1_v[i] + 2.0 * k2_v[i] + 2.0 * k3_v[i] + k4_v[i]) / 6.0; + } +} + +VelocityVerletSolver::VelocityVerletSolver(PenningTrap& trap_in, double dt_in) + : Solver(trap_in, dt_in) { + // Initialize accelerations + int n = trap.num_particles(); + accelerations.resize(n); + for (int i = 0; i < n; i++) { + Particle& p = trap.get_particle(i); + arma::vec F = trap.total_force(i, 0.0); + accelerations[i] = F / p.mass; + } +} + +// Perform a single Velocity Verlet integration step +void VelocityVerletSolver::step() { + int n = trap.num_particles(); + for (int i = 0; i < n; i++) { + Particle& p = trap.get_particle(i); + p.position += p.velocity * dt + 0.5 * accelerations[i] * dt * dt; + arma::vec new_acceleration = trap.total_force(i, steps * dt) / p.mass; + p.velocity += 0.5 * (accelerations[i] + new_acceleration) * dt; + accelerations[i] = new_acceleration; + } +} + + +// Initialize Analytical solver +AnalyticalSolver::AnalyticalSolver(PenningTrap& trap_in, double dt_in) + : Solver(trap_in, dt_in) { + // Check for validity of analytical solution + if (trap_in.get_coulomb_interaction()) { + throw std::runtime_error("Analytical solution not valid with Coulomb interactions enabled."); + } + if (trap_in.get_modified_fields()) { + throw std::runtime_error("Analytical solution not valid with modified external fields."); + } + for (int i = 0; i < trap_in.num_particles(); i++) { + if (trap_in.get_particle(i).velocity(0) != 0 || trap_in.get_particle(i).velocity(2) != 0) { + throw std::runtime_error("Analytical solution not valid for non-zero initial velocity in x or z direction."); + } + if (trap_in.get_particle(i).position(1) != 0) { + throw std::runtime_error("Analytical solution not valid for non-zero initial position in y direction."); + } + } + + trap = trap_in; + dt = dt_in * constants::second; + steps = 0; + int n = trap.num_particles(); + w_0.resize(n); + w_z.resize(n); + w_plus.resize(n); + w_minus.resize(n); + A_plus.resize(n); + A_minus.resize(n); + z_0.resize(n); + + for (int i = 0; i < n; i++) { + Particle& p = trap.get_particle(i); + double q = p.charge; + double m = p.mass; + double B0 = trap.B0; + double V0 = trap.V0; + double d = trap.d; + + w_0[i] = q * B0 / m; + w_z[i] = std::sqrt(2 * q * V0 / (m * d * d)); + double discriminant = w_0[i] * w_0[i] - 2 * w_z[i] * w_z[i]; + if (discriminant < 0) { + throw std::runtime_error("Analytical solution not valid: discriminant is negative."); + } + w_plus[i] = (w_0[i] + std::sqrt(discriminant)) / 2.0; + w_minus[i] = (w_0[i] - std::sqrt(discriminant)) / 2.0; + + arma::vec r = p.position; + arma::vec v = p.velocity; + + // Initial conditions + double x0 = r(0); + double y0 = r(1); + double z0 = r(2); + double vx0 = v(0); + double vy0 = v(1); + double vz0 = v(2); + + // Calculate amplitudes + A_plus[i] = (vy0 + w_minus[i] * x0) / (w_plus[i] - w_minus[i]); + A_minus[i] = -(vy0 + w_plus[i] * x0) / (w_plus[i] - w_minus[i]); + z_0[i] = z0; + } +} + +// Perform a single Analytical "integration" step +void AnalyticalSolver::step() { + int n = trap.num_particles(); + double t = steps * dt; + for (int i = 0; i < n; i++) { + Particle& p = trap.get_particle(i); + double x = -A_plus[i] * std::cos(w_plus[i] * t) - A_minus[i] * std::cos(w_minus[i] * t); + double y = A_plus[i] * std::sin(w_plus[i] * t) + A_minus[i] * std::sin(w_minus[i] * t); + double z = z_0[i] * std::cos(w_z[i] * t); + + double vx = +A_plus[i] * w_plus[i] * std::sin(w_plus[i] * t) + A_minus[i] * w_minus[i] * std::sin(w_minus[i] * t); + double vy = A_plus[i] * w_plus[i] * std::cos(w_plus[i] * t) + A_minus[i] * w_minus[i] * std::cos(w_minus[i] * t); + double vz = -z_0[i] * w_z[i] * std::sin(w_z[i] * t); + + p.position = arma::vec({x, y, z}); + p.velocity = arma::vec({vx, vy, vz}); + } +} + + +// Convert vector of matrices to a 3D cube for easier data handling +arma::cube to_cube(const std::vector& vec) { + if (vec.empty()) { + return arma::cube(); + } + int rows = vec[0].n_rows; + int cols = vec[0].n_cols; + int slices = vec.size(); + arma::cube cube_data(rows, cols, slices); + for (size_t i = 0; i < vec.size(); ++i) { + cube_data.slice(i) = vec[i]; + } + return cube_data; +} + +// Save cube data to a file +void save_to_file(const std::string& filename, const arma::cube& data) { + data.save(filename, arma::raw_ascii); +} + +void save_to_file(const std::string& filename, const arma::vec& data) { + arma::mat mat_data = arma::conv_to::from(data); + mat_data.save(filename, arma::raw_ascii); +} \ No newline at end of file diff --git a/src/project3/two_particles.cpp b/src/project3/two_particles.cpp new file mode 100644 index 0000000..f247d2d --- /dev/null +++ b/src/project3/two_particles.cpp @@ -0,0 +1,126 @@ +#include +#include "classes.hpp" +#include "solvers.hpp" +#include "constants.hpp" +#include "argparse/argparse.hpp" + +#include + +using namespace std; + +int main(int argc, char* argv[]) { + argparse::ArgumentParser program("two_particles"); + program.add_argument("-N", "--steps") + .help("Number of time steps") + .default_value(4000) + .scan<'i', int>(); + program.add_argument("-n", "--num_particles") + .help("Number of particles (either 1 or 2)") + .default_value(2) + .scan<'i', int>(); + program.add_argument("-i", "--disable-interactions") + .help("Disable Coulomb interactions") + .default_value(true) + .implicit_value(false); + program.add_argument("-E", "--euler") + .help("Use Euler solver instead of RK4") + .flag(); + program.add_argument("-A", "--analytical") + .help("Use Analytical solver instead of RK4 (only valid for particles without interactions and specific initial conditions)") + .flag(); + program.add_argument("-V", "--velocity-verlet") + .help("Use Velocity Verlet solver instead of RK4") + .flag(); + + try { + program.parse_args(argc, argv); + } catch (const runtime_error& err) { + cerr << err.what() << endl; + cerr << program.help().str() << endl; + exit(1); + } + + double N = program.get("--steps"); + int n_particles = program.get("--num_particles"); + bool interactions = program.get("--disable-interactions"); + bool use_euler = program.get("--euler"); + bool use_analytical = program.get("--analytical"); + bool use_velocity_verlet = program.get("--velocity-verlet"); + + double B_0 = 1.0; // Tesla + double V_0 = 0.025; // Volt + double d = 500e-6; // meter + + PenningTrap trap(B_0, V_0, d); + if (interactions) { + trap.enable_interactions(); + } + else { + trap.disable_interactions(); + } + cout << trap << endl; + + + Particle p1(arma::vec({20e-6, 0, 20e-6}), arma::vec({0, 25, 0}), 40.078 * constants::amu, 1 * constants::elementary_charge); + Particle p2(arma::vec({25e-6, 25e-6, 0}), arma::vec({0, 40, 5}), 40.078 * constants::amu, 1 * constants::elementary_charge); + + trap.add_particle(p1); + if (n_particles > 1) { + trap.add_particle(p2); + } + + double dt = 50e-6 / N; // seconds + + if (use_euler) { + EulerSolver solver(trap, dt); + solver.simulate(N); + solver.save("results/two_particles_" + std::to_string(n_particles) + (interactions ? "_with_interactions_" : "_no_interactions_") + std::to_string(N) + "_steps_euler"); + vector positions = solver.get_positions(); + cout << "Final position of particle 1: " << endl; + cout << positions.back().col(0) << endl; + if (n_particles > 1) { + cout << "Final position of particle 2: " << endl; + cout << positions.back().col(1) << endl; + } + return 0; + } + + if (use_analytical) { + AnalyticalSolver solver(trap, dt); + solver.simulate(N); + solver.save("results/two_particles_" + std::to_string(n_particles) + (interactions ? "_with_interactions_" : "_no_interactions_") + std::to_string(N) + "_steps_analytical"); + vector positions = solver.get_positions(); + cout << "Final position of particle 1: " << endl; + cout << positions.back().col(0) << endl; + if (n_particles > 1) { + cout << "Final position of particle 2: " << endl; + cout << positions.back().col(1) << endl; + } + return 0; + } + + if (use_velocity_verlet) { + VelocityVerletSolver solver(trap, dt); + solver.simulate(N); + solver.save("results/two_particles_" + std::to_string(n_particles) + (interactions ? "_with_interactions_" : "_no_interactions_") + std::to_string(N) + "_steps_velocity_verlet"); + vector positions = solver.get_positions(); + cout << "Final position of particle 1: " << endl; + cout << positions.back().col(0) << endl; + if (n_particles > 1) { + cout << "Final position of particle 2: " << endl; + cout << positions.back().col(1) << endl; + } + return 0; + } + + RK4Solver solver(trap, dt); + solver.simulate(N); + solver.save("results/two_particles_" + std::to_string(n_particles) + (interactions ? "_with_interactions_" : "_no_interactions_") + std::to_string(N) + "_steps_rk4"); + vector positions = solver.get_positions(); + cout << "Final position of particle 1: " << endl; + cout << positions.back().col(0) << endl; + if (n_particles > 1) { + cout << "Final position of particle 2: " << endl; + cout << positions.back().col(1) << endl; + } +} \ No newline at end of file diff --git a/uv.lock b/uv.lock index ce9600a..c211472 100644 --- a/uv.lock +++ b/uv.lock @@ -42,6 +42,18 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/7c/fc/6a8cb64e5f0324877d503c854da15d76c1e50eb722e320b15345c4d0c6de/cffi-1.17.1-cp313-cp313-win_amd64.whl", hash = "sha256:f6a16c31041f09ead72d69f583767292f750d24913dadacf5756b966aacb3f1a", size = 182009, upload-time = "2024-09-04T20:44:45.309Z" }, ] +[[package]] +name = "click" +version = "8.3.0" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "colorama", marker = "sys_platform == 'win32'" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/46/61/de6cd827efad202d7057d93e0fed9294b96952e188f7384832791c7b2254/click-8.3.0.tar.gz", hash = "sha256:e7b8232224eba16f4ebe410c25ced9f7875cb5f3263ffc93cc3e8da705e229c4", size = 276943, upload-time = "2025-09-18T17:32:23.696Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/db/d3/9dcc0f5797f070ec8edf30fbadfb200e71d9db6b84d211e3b2085a7589a0/click-8.3.0-py3-none-any.whl", hash = "sha256:9b9f285302c6e3064f4330c05f05b81945b2a39544279343e6e7c5f27a9baddc", size = 107295, upload-time = "2025-09-18T17:32:22.42Z" }, +] + [[package]] name = "colorama" version = "0.4.6" @@ -198,6 +210,7 @@ dependencies = [ { name = "numpy" }, { name = "pandas" }, { name = "scipy" }, + { name = "typer" }, ] [package.metadata] @@ -207,6 +220,7 @@ requires-dist = [ { name = "numpy", specifier = ">=2.3.2" }, { name = "pandas", specifier = ">=2.3.2" }, { name = "scipy", specifier = ">=1.16.1" }, + { name = "typer", specifier = ">=0.19.2" }, ] [[package]] @@ -367,6 +381,18 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/80/be/3578e8afd18c88cdf9cb4cffde75a96d2be38c5a903f1ed0ceec061bd09e/kiwisolver-1.4.9-cp314-cp314t-win_arm64.whl", hash = "sha256:4a48a2ce79d65d363597ef7b567ce3d14d68783d2b2263d98db3d9477805ba32", size = 70260, upload-time = "2025-08-10T21:27:36.606Z" }, ] +[[package]] +name = "markdown-it-py" +version = "4.0.0" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "mdurl" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/5b/f5/4ec618ed16cc4f8fb3b701563655a69816155e79e24a17b651541804721d/markdown_it_py-4.0.0.tar.gz", hash = "sha256:cb0a2b4aa34f932c007117b194e945bd74e0ec24133ceb5bac59009cda1cb9f3", size = 73070, upload-time = "2025-08-11T12:57:52.854Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/94/54/e7d793b573f298e1c9013b8c4dade17d481164aa517d1d7148619c2cedbf/markdown_it_py-4.0.0-py3-none-any.whl", hash = "sha256:87327c59b172c5011896038353a81343b6754500a08cd7a4973bb48c6d578147", size = 87321, upload-time = "2025-08-11T12:57:51.923Z" }, +] + [[package]] name = "matplotlib" version = "3.10.5" @@ -426,6 +452,15 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/8f/8e/9ad090d3553c280a8060fbf6e24dc1c0c29704ee7d1c372f0c174aa59285/matplotlib_inline-0.1.7-py3-none-any.whl", hash = "sha256:df192d39a4ff8f21b1895d72e6a13f5fcc5099f00fa84384e0ea28c2cc0653ca", size = 9899, upload-time = "2024-04-15T13:44:43.265Z" }, ] +[[package]] +name = "mdurl" +version = "0.1.2" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/d6/54/cfe61301667036ec958cb99bd3efefba235e65cdeb9c84d24a8293ba1d90/mdurl-0.1.2.tar.gz", hash = "sha256:bb413d29f5eea38f31dd4754dd7377d4465116fb207585f97bf925588687c1ba", size = 8729, upload-time = "2022-08-14T12:40:10.846Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/b3/38/89ba8ad64ae25be8de66a6d463314cf1eb366222074cfda9ee839c56a4b4/mdurl-0.1.2-py3-none-any.whl", hash = "sha256:84008a41e51615a49fc9966191ff91509e3c40b939176e643fd50a5c2196b8f8", size = 9979, upload-time = "2022-08-14T12:40:09.779Z" }, +] + [[package]] name = "nest-asyncio" version = "1.6.0" @@ -757,6 +792,19 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/69/a4/7ee652ea1c77d872f5d99ed937fa8bbd1f6f4b7a39a6d3a0076c286e0c3e/pyzmq-27.0.2-cp314-cp314t-win_arm64.whl", hash = "sha256:4108785f2e5ac865d06f678a07a1901e3465611356df21a545eeea8b45f56265", size = 574901, upload-time = "2025-08-21T04:22:17.423Z" }, ] +[[package]] +name = "rich" +version = "14.1.0" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "markdown-it-py" }, + { name = "pygments" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/fe/75/af448d8e52bf1d8fa6a9d089ca6c07ff4453d86c65c145d0a300bb073b9b/rich-14.1.0.tar.gz", hash = "sha256:e497a48b844b0320d45007cdebfeaeed8db2a4f4bcf49f15e455cfc4af11eaa8", size = 224441, upload-time = "2025-07-25T07:32:58.125Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/e3/30/3c4d035596d3cf444529e0b2953ad0466f6049528a879d27534700580395/rich-14.1.0-py3-none-any.whl", hash = "sha256:536f5f1785986d6dbdea3c75205c473f970777b4a0d6c6dd1b696aa05a3fa04f", size = 243368, upload-time = "2025-07-25T07:32:56.73Z" }, +] + [[package]] name = "scipy" version = "1.16.1" @@ -804,6 +852,15 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/6e/6c/a76329897a7cae4937d403e623aa6aaea616a0bb5b36588f0b9d1c9a3739/scipy-1.16.1-cp314-cp314t-win_amd64.whl", hash = "sha256:c0c804d60492a0aad7f5b2bb1862f4548b990049e27e828391ff2bf6f7199998", size = 39427705, upload-time = "2025-07-27T16:31:53.96Z" }, ] +[[package]] +name = "shellingham" +version = "1.5.4" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/58/15/8b3609fd3830ef7b27b655beb4b4e9c62313a4e8da8c676e142cc210d58e/shellingham-1.5.4.tar.gz", hash = "sha256:8dbca0739d487e5bd35ab3ca4b36e11c4078f3a234bfce294b0a0291363404de", size = 10310, upload-time = "2023-10-24T04:13:40.426Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/e0/f9/0595336914c5619e5f28a1fb793285925a8cd4b432c9da0a987836c7f822/shellingham-1.5.4-py2.py3-none-any.whl", hash = "sha256:7ecfff8f2fd72616f7481040475a65b2bf8af90a56c89140852d1120324e8686", size = 9755, upload-time = "2023-10-24T04:13:38.866Z" }, +] + [[package]] name = "six" version = "1.17.0" @@ -855,6 +912,30 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/00/c0/8f5d070730d7836adc9c9b6408dec68c6ced86b304a9b26a14df072a6e8c/traitlets-5.14.3-py3-none-any.whl", hash = "sha256:b74e89e397b1ed28cc831db7aea759ba6640cb3de13090ca145426688ff1ac4f", size = 85359, upload-time = "2024-04-19T11:11:46.763Z" }, ] +[[package]] +name = "typer" +version = "0.19.2" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "click" }, + { name = "rich" }, + { name = "shellingham" }, + { name = "typing-extensions" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/21/ca/950278884e2ca20547ff3eb109478c6baf6b8cf219318e6bc4f666fad8e8/typer-0.19.2.tar.gz", hash = "sha256:9ad824308ded0ad06cc716434705f691d4ee0bfd0fb081839d2e426860e7fdca", size = 104755, upload-time = "2025-09-23T09:47:48.256Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/00/22/35617eee79080a5d071d0f14ad698d325ee6b3bf824fc0467c03b30e7fa8/typer-0.19.2-py3-none-any.whl", hash = "sha256:755e7e19670ffad8283db353267cb81ef252f595aa6834a0d1ca9312d9326cb9", size = 46748, upload-time = "2025-09-23T09:47:46.777Z" }, +] + +[[package]] +name = "typing-extensions" +version = "4.15.0" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/72/94/1a15dd82efb362ac84269196e94cf00f187f7ed21c242792a923cdb1c61f/typing_extensions-4.15.0.tar.gz", hash = "sha256:0cea48d173cc12fa28ecabc3b837ea3cf6f38c6d1136f85cbaaf598984861466", size = 109391, upload-time = "2025-08-25T13:49:26.313Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/18/67/36e9267722cc04a6b9f15c7f3441c2363321a3ea07da7ae0c0707beb2a9c/typing_extensions-4.15.0-py3-none-any.whl", hash = "sha256:f0fa19c6845758ab08074a0cfa8b7aecb71c999ca73d62883bc25cc018c4e548", size = 44614, upload-time = "2025-08-25T13:49:24.86Z" }, +] + [[package]] name = "tzdata" version = "2025.2"