From ff8454900b142d24eea96abb7c671feae312e07c Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Tue, 12 May 2020 00:47:42 -0400
Subject: Two changes to MCMC library.   - Added orthogonalize command that
 resets matrix to a nearby     orthogonal one.   - run command now returns
 proportion of acceptances.

---
 hadamard_mcmc.hpp | 39 +++++++++++++++++++++++++++++++++++++--
 1 file changed, 37 insertions(+), 2 deletions(-)

diff --git a/hadamard_mcmc.hpp b/hadamard_mcmc.hpp
index bb948de..41d7a48 100644
--- a/hadamard_mcmc.hpp
+++ b/hadamard_mcmc.hpp
@@ -40,6 +40,34 @@ public:
 
     return E;
   }
+
+  void orthogonalize() {
+    for (unsigned i = 0; i < d; i++) {
+      for (unsigned j = 0; j < i; j++) {
+        double ij = 0;
+        for (unsigned k = 0; k < d; k++) {
+          ij += (*this)(i, k) * (*this)(j, k);
+        }
+
+        ij /= d;
+
+        for (unsigned k = 0; k < d; k++) {
+          (*this)(i, k) -= ij * (*this)(j, k);
+        }
+      }
+
+      double ii = 0;
+      for (unsigned j = 0; j < d; j++) {
+        ii += pow((*this)(i, j), 2);
+      }
+
+      ii = sqrt(ii / d);
+
+      for (unsigned j = 0; j < d; j++) {
+        (*this)(i, j) = (*this)(i, j) / ii;
+      }
+    }
+  }
 };
 
 class Givens {
@@ -187,13 +215,20 @@ public:
     }
   }
 
-  void run(unsigned N, bool dry = false) {
+  double run(unsigned N, bool dry = false, unsigned nOrth = 1000) {
+    double totalRatio = 0;
+
     for (unsigned i = 0; i < N; i++) {
-      this->sweep(dry);
+      totalRatio += this->sweep(dry);
       if (!dry) {
         A.after_sweep(θ0, E, M);
       }
+      if (i % nOrth == 0) {
+        M.orthogonalize();
+      }
     }
+
+    return totalRatio / N;
   }
 };
 
-- 
cgit v1.2.3-70-g09d2