diff options
| -rw-r--r-- | langevin.cpp | 19 |
1 files 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; |