blob: 3683c10e39ceff8fe1728311069096283a2e8077 (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
|
#include "mp.h"
#include <algorithm>
#include <vector>
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 = (arma::uword) Y.n_rows;
V i(n), j(n);
for (arma::uword k = 0; k < n; k++)
i[k] = j[k] = k;
double prev_delta_sum = 1. / 0.;
for (size_t iter = 0; iter < maxIter; iter++) {
double delta_sum = 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), mp::EPSILON);
double delta = (D(*a, *b) - d2) / fraction;
delta_sum += fabs(delta);
Y.row(*b) += delta * (direction / d2);
}
}
if (fabs(prev_delta_sum - delta_sum) < tol)
break;
prev_delta_sum = delta_sum;
}
return Y;
}
|