Add Random walk example, Improve vec class

Makefile

@@ -2,7 +2,7 @@ default_target: all
 .PHONY: all example test_core test_algebra test_polynomial test_optimization test_autodiff \
 		test_calculus test_statistics automatic_differentiation hamiltonian error_propagation \
 		statistics sampling_distributions multivariate_minimization benchmark \
-		clean attractor montecarlo_comparison histogram
+		clean attractor montecarlo_comparison histogram random_walk
 
 all: examples test
 
@@ -116,9 +116,13 @@ histogram:
 	@echo Compiling \"Histogram\" example...
 	@g++ examples/histogram.cpp ${CXXFLAGS} -o ./examples/histogram
 
+random_walk:
+	@echo Compiling \"Random walk\" example...
+	@g++ examples/random_walk.cpp ${CXXFLAGS} -o ./examples/random_walk
+
 examples: example automatic_differentiation hamiltonian error_propagation \
 		  statistics sampling_distributions montecarlo_comparison multivariate_minimization \
-		  logfit attractor histogram
+		  logfit attractor histogram random_walk
 
 
 # Benchmarks

examples/random_walk.cpp

@@ -0,0 +1,73 @@
+
+///
+/// @file random_walk.cpp A 3D random walk generator
+/// This example may be compiled using 'make random_walk'
+///
+
+#include "theoretica.h"
+using namespace th;
+
+#include <fstream>
+#include <ctime>
+#include <algorithm>
+
+
+// You can plot the data file using gnuplot:
+// splot "random_walk.dat" with lines title "3D Random Walk"
+
+
+int main(int argc, char* argv[]) {
+
+    // Number of steps
+    const unsigned int N = 10000;
+
+    // Output file
+    std::ofstream file("random_walk.dat");
+
+    // Random number generator
+    PRNG g = PRNG::xoshiro(time(nullptr));
+
+    // Trajectory
+    std::vector<vec3> pos(N);
+
+    // Start from the origin
+    pos[0] = {0, 0, 0};
+    file << pos[0].to_string(" ", false) << std::endl;
+
+    for (unsigned int i = 1; i < N; ++i) {
+
+        // Generate a random vector with fixed length
+        // You can change the generation of r to test
+        // different distributions (e.g. rand_gaussian(0, 1, g))
+        const real r = 1;
+        const real theta = rand_uniform(0, PI, g);
+        const real phi = rand_uniform(0, TAU, g);
+
+        // Update the trajectory
+        pos[i] = pos[i - 1] + vec3({
+            r * sin(theta) * cos(phi),
+            r * sin(theta) * sin(phi),
+            r * cos(theta)
+        });
+
+        file << pos[i].to_string(" ", false) << std::endl;
+    }
+
+
+    // Compute the mean direction of the steps
+    vec3 mean_dir = {0, 0, 0};
+    for (unsigned int i = 1; i < N; ++i)
+        mean_dir += (pos[i] - pos[i - 1]) / N;
+
+    std::cout << "Mean Direction: " << mean_dir << std::endl;
+
+
+    // Compute the RMS displacement
+    real sqr_disp = 0;
+    for (unsigned int i = 0; i < N; ++i)
+        sqr_disp += pos[i] * pos[i] / N;
+
+    std::cout << "RMS Displacement: " << sqrt(sqr_disp) << std::endl;
+
+    return 0;
+}

src/algebra/vec.h

@@ -311,17 +311,23 @@ namespace theoretica {
 #ifndef THEORETICA_NO_PRINT
 
 		/// Convert the vector to string representation
-		inline std::string to_string(const std::string& separator = ", ") const {
+		inline std::string to_string(
+			const std::string& separator = ", ",
+			bool parenthesis = true) const {
 
 			std::stringstream res;
 
-			res << "(";
+			if(parenthesis)
+				res << "(";
+
 			for (unsigned int i = 0; i < N; ++i) {
 				res << data[i];
 				if(i != N - 1)
 					res << separator;
 			}
-			res << ")";
+
+			if(parenthesis)
+				res << ")";
 
 			return res.str();
 		}
@@ -653,17 +659,23 @@ namespace theoretica {
 #ifndef THEORETICA_NO_PRINT
 
 		/// Convert the vector to string representation
-		inline std::string to_string(const std::string& separator = ", ") const {
+		inline std::string to_string(
+			const std::string& separator = ", ",
+			bool parenthesis = true) const {
 
 			std::stringstream res;
 
-			res << "(";
+			if(parenthesis)
+				res << "(";
+
 			for (unsigned int i = 0; i < size(); ++i) {
 				res << data[i];
 				if(i != size() - 1)
 					res << separator;
 			}
-			res << ")";
+
+			if(parenthesis)
+				res << ")";
 
 			return res.str();
 		}