From 07c5fee22dd83f791b97567df301cacc3f147d23 Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Thu, 26 Dec 2024 11:34:59 +0100
Subject: Another small refactor.

---
 walk.cpp | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/walk.cpp b/walk.cpp
index f064bb5..45b3601 100644
--- a/walk.cpp
+++ b/walk.cpp
@@ -24,6 +24,10 @@ Vector randomVector(unsigned N, Rng& r, Real σ = 1) {
   return v;
 }
 
+Vector projectionOfOn(const Vector& v, const Vector& u) {
+  return (v.dot(u) / u.squaredNorm()) * u;
+}
+
 Real wignerCDF(Real λ) {
   return 0.5 + (λ * sqrt(4 - pow(λ, 2)) / 4 + atan(λ / sqrt(4 - pow(λ, 2)))) / M_PI;
 }
@@ -63,7 +67,7 @@ public:
 
   Vector ∇H(const Vector& x) const {
     Vector ∂H = J.cwiseProduct(x);
-    return ∂H - (∂H.dot(x) / x.squaredNorm()) * x;
+    return ∂H - projectionOfOn(∂H, x);
   }
 };
 
@@ -98,9 +102,8 @@ Vector gradientDescent(const QuadraticModel& M, const Vector& x₀, Real E, Real
 
 Vector randomStep(const QuadraticModel& M, const Vector& x₀, Real E, Rng& r, Real Δt = 1e-4) {
   Vector η = randomVector(M.N, r, sqrt(2 * Δt));
-  η -= (η.dot(x₀) / x₀.squaredNorm()) * x₀;
-  Vector ∇H₀ = M.∇H(x₀);
-  η -= (η.dot(∇H₀) / ∇H₀.squaredNorm()) * ∇H₀;
+  η -= projectionOfOn(η, x₀);
+  η -= projectionOfOn(η, M.∇H(x₀));
   return gradientDescent(M, normalizeVector(x₀ + η), E);
 }
 
-- 
cgit v1.2.3-70-g09d2