From 8444cda9634e53aa1987a0f5039f73c815b67d3b Mon Sep 17 00:00:00 2001
From: Samuel Fadel <samuelfadel@gmail.com>
Date: Fri, 29 Jan 2016 10:12:33 +0100
Subject: Renamed ForceScheme source file to lowercase.

---
 forcescheme.cpp | 48 ++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 48 insertions(+)
 create mode 100644 forcescheme.cpp

(limited to 'forcescheme.cpp')

diff --git a/forcescheme.cpp b/forcescheme.cpp
new file mode 100644
index 0000000..22ab293
--- /dev/null
+++ b/forcescheme.cpp
@@ -0,0 +1,48 @@
+#include "mp.h"
+
+#include <algorithm>
+
+static const double EPSILON = 1e-3;
+
+typedef arma::uvec V;
+
+arma::mat mp::forceScheme(const arma::mat &D,
+        arma::mat &Y,
+        size_t maxIter,
+        double tol,
+        double fraction)
+{
+    arma::uword n = Y.n_rows;
+    V i(n), j(n);
+    for (arma::uword k = 0; k < n; k++) {
+        i[k] = j[k] = k;
+    }
+
+    double prevDeltaSum = 1. / 0.;
+    for (size_t iter = 0; iter < maxIter; iter++) {
+        double deltaSum = 0;
+
+        arma::shuffle(i);
+        for (V::iterator a = i.begin(); a != i.end(); a++) {
+            arma::shuffle(j);
+            for (V::iterator b = j.begin(); b != j.end(); b++) {
+                if (*a == *b) {
+                    continue;
+                }
+
+                arma::rowvec direction(Y.row(*b) - Y.row(*a));
+                double d2 = std::max(arma::norm(direction, 2), EPSILON);
+                double delta = (D(*a, *b) - d2) / fraction;
+                deltaSum += fabs(delta);
+                Y.row(*b) += delta * (direction / d2);
+            }
+        }
+
+        if (fabs(prevDeltaSum - deltaSum) < tol) {
+            break;
+        }
+        prevDeltaSum = deltaSum;
+    }
+
+    return Y;
+}
-- 
cgit v1.2.3