From 7c3e71970b6f2d48530bc2ab4fc0f2932522b98b Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Tue, 5 Jan 2021 12:18:56 +0100
Subject: New function to create random vectors.

---
 langevin.cpp | 19 +++++++------------
 1 file changed, 7 insertions(+), 12 deletions(-)

diff --git a/langevin.cpp b/langevin.cpp
index e4b6c8d..9da4ef6 100644
--- a/langevin.cpp
+++ b/langevin.cpp
@@ -12,16 +12,14 @@
 
 using Rng = randutils::random_generator<pcg32>;
 
-Vector initializeVector(unsigned N, double a, Rng& r) {
+template <class Distribution, class Generator>
+Vector randomVector(unsigned N, Distribution d, Generator& r) {
   Vector z(N);
-  complex_normal_distribution<> dist(0, a, 0);
 
   for (unsigned i = 0; i < N; i++) {
-    z(i) = dist(r.engine());
+    z(i) = d(r);
   }
 
-  z *= sqrt(N) / sqrt(z.dot(z)); // Normalize.
-
   return z;
 }
 
@@ -38,10 +36,7 @@ Vector langevin(const Tensor& J, const Vector& z0, double T, double γ0,
 
   while (!quit(W, steps)) {
     double γ = pow(r.variate<double, std::normal_distribution>(0, γ0), 2);
-    Vector η(z.size());
-    for (unsigned i = 0; i < z.size(); i++) {
-      η(i) = d(r.engine());
-    }
+    Vector η = randomVector(z.size(), d, r.engine());
 
     Vector zNext = z - γ * dW + η;
     zNext *= sqrt(zNext.size()) / sqrt(zNext.dot(zNext));
@@ -105,10 +100,10 @@ int main(int argc, char* argv[]) {
   double σ = sqrt(factorial(p) / (2.0 * pow(N, p - 1)));
 
   Rng r;
-  complex_normal_distribution<> d(0, σ, κ);
 
-  Tensor J = generateCouplings<Scalar, PSPIN_P>(N, d, r.engine());
-  Vector z = initializeVector(N, 100, r);
+  Tensor J = generateCouplings<Scalar, PSPIN_P>(N, complex_normal_distribution<>(0, σ, κ), r.engine());
+  Vector z = randomVector(N, complex_normal_distribution<>(0, 1, 0), r.engine());
+  z *= sqrt(N) / sqrt(z.dot(z)); // Normalize.
 
   std::function<bool(double, unsigned)> findSaddle = [δ](double W, unsigned) {
     std::cout << W << std::endl;
-- 
cgit v1.2.3-70-g09d2