main(): fixed code for init'ing CP projection.

Minor code improvements in forceScheme() and dist().

3 files changed, 23 insertions, 23 deletions

diff --git a/dist.cpp b/dist.cpp
index f307564..cd85b6f 100644
--- a/dist.cpp
+++ b/dist.cpp
@@ -10,11 +10,12 @@ double mp::euclidean(const arma::rowvec &x1, const arma::rowvec &x2)
 arma::mat mp::dist(const arma::mat &X, mp::DistFunc dfunc)
 {
     int n = uintToInt<arma::uword, int>(X.n_rows);
-    arma::mat D(n, n, arma::fill::zeros);
+    arma::mat D(n, n);
 
     #pragma omp parallel for shared(X, D, n)
     for (int i = 0; i < n; i++) {
-        for (arma::uword j = 0; j < i; j++) {
+        D(i, i) = 0;
+        for (int j = 0; j < i; j++) {
             D(i, j) = dfunc(X.row(i), X.row(j));
             D(j, i) = D(i, j);
         }
diff --git a/forcescheme.cpp b/forcescheme.cpp
index 9e3feaa..44c3acc 100644
--- a/forcescheme.cpp
+++ b/forcescheme.cpp
@@ -5,8 +5,6 @@
 
 static const double EPSILON = 1e-3;
 
-typedef arma::uvec V;
-
 arma::mat mp::forceScheme(const arma::mat &D,
         arma::mat &Y,
         size_t maxIter,
@@ -14,28 +12,28 @@ arma::mat mp::forceScheme(const arma::mat &D,
         double fraction)
 {
     arma::uword n = Y.n_rows;
-    V i(n), j(n);
+    arma::uvec indices1(n), indices2(n);
     for (arma::uword k = 0; k < n; k++) {
-        i[k] = j[k] = k;
+        indices1[k] = indices2[k] = k;
     }
 
     double prevDeltaSum = std::numeric_limits<double>::infinity();
     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) {
+        arma::shuffle(indices1);
+        for (auto i: indices1) {
+            arma::shuffle(indices2);
+            for (auto j: indices2) {
+                if (i == j) {
                     continue;
                 }
 
-                arma::rowvec direction(Y.row(*b) - Y.row(*a));
+                arma::rowvec direction(Y.row(j) - Y.row(i));
                 double d2 = std::max(arma::norm(direction, 2), EPSILON);
-                double delta = (D(*a, *b) - d2) / fraction;
+                double delta = (D(i, j) - d2) / fraction;
                 deltaSum += fabs(delta);
-                Y.row(*b) += delta * (direction / d2);
+                Y.row(j) += delta * (direction / d2);
             }
         }
 
diff --git a/main.cpp b/main.cpp
index 7a301a7..e778dca 100644
--- a/main.cpp
+++ b/main.cpp
@@ -23,7 +23,7 @@
 #include "selectionhandler.h"
 #include "brushinghandler.h"
 
-static const int RNG_SEED = 1;
+static const int RNG_SEED = 96123;
 
 static QObject *mainProvider(QQmlEngine *engine, QJSEngine *scriptEngine)
 {
@@ -36,12 +36,11 @@ static QObject *mainProvider(QQmlEngine *engine, QJSEngine *scriptEngine)
 arma::uvec extractCPs(const arma::mat &X)
 {
     int numCPs = (int) (3 * sqrt(X.n_rows));
-    arma::uvec cpIndices(numCPs);
-    std::iota(cpIndices.begin(), cpIndices.end(), 0);
-    std::shuffle(cpIndices.begin(),
-                 cpIndices.end(),
-                 std::default_random_engine(RNG_SEED));
-    return cpIndices;
+    arma::uvec indices(X.n_rows);
+    std::iota(indices.begin(), indices.end(), 0);
+    arma::shuffle(indices);
+
+    return indices.subvec(0, numCPs-1);
 }
 
 int main(int argc, char **argv)
@@ -80,9 +79,10 @@ int main(int argc, char **argv)
         return 1;
     }
 
-    arma::mat X = m->X();
-    arma::vec labels = m->labels();
+    const arma::mat &X = m->X();
+    //arma::vec labels = m->labels();
 
+    arma::arma_rng::set_seed(RNG_SEED);
     arma::uvec cpIndices;
     arma::mat Ys;
 
@@ -112,6 +112,7 @@ int main(int argc, char **argv)
     } else {
         std::cerr << "No CP file, generating initial projection...\n";
         Ys.set_size(cpIndices.n_elem, 2);
+        Ys.randu();
         mp::forceScheme(mp::dist(X.rows(cpIndices)), Ys);
     }