aboutsummaryrefslogtreecommitdiff
path: root/forcescheme.cpp
diff options
context:
space:
mode:
authorSamuel Fadel <samuelfadel@gmail.com>2016-01-29 10:12:33 +0100
committerSamuel Fadel <samuelfadel@gmail.com>2016-01-29 10:12:56 +0100
commit8444cda9634e53aa1987a0f5039f73c815b67d3b (patch)
tree6b2f51d11cbdf8bd7a97f1db1fd7487a4cc5c4e0 /forcescheme.cpp
parent990a1eb8aa57aecddf0b8d0d8c230875a7421a10 (diff)
Renamed ForceScheme source file to lowercase.
Diffstat (limited to 'forcescheme.cpp')
-rw-r--r--forcescheme.cpp48
1 files changed, 48 insertions, 0 deletions
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;
+}