Initial commit. ForceScheme seems bugged.

12 files changed, 466 insertions, 0 deletions

diff --git a/dist.cpp b/dist.cpp
new file mode 100644
index 0000000..aaa2167
--- /dev/null
+++ b/dist.cpp
@@ -0,0 +1,20 @@
+#include "mp.h"
+
+double mp::euclidean(const arma::rowvec &x1, const arma::rowvec &x2)
+{
+    return arma::norm(x1 - x2, 2);
+}
+
+arma::mat mp::dist(const arma::mat &X, double (*distCalc)(const arma::rowvec &, const arma::rowvec &))
+{
+    arma::uword n = X.n_rows;
+    arma::mat D(n, n, arma::fill::zeros);
+
+    for (arma::uword i = 0; i < n; i++)
+        for (arma::uword j = 0; j < i; j++) {
+            D(i, j) = distCalc(X.row(i), X.row(j));
+            D(j, i) = D(i, j);
+        }
+
+    return D;
+}
diff --git a/forceScheme.cpp b/forceScheme.cpp
new file mode 100644
index 0000000..3683c10
--- /dev/null
+++ b/forceScheme.cpp
@@ -0,0 +1,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;
+}
diff --git a/glyph.cpp b/glyph.cpp
new file mode 100644
index 0000000..f143ee8
--- /dev/null
+++ b/glyph.cpp
@@ -0,0 +1,103 @@
+#include "glyph.h"
+
+#include <cmath>
+#include <QSGGeometry>
+#include <QSGGeometryNode>
+#include <QSGMaterial>
+#include <QSGFlatColorMaterial>
+
+const float PI = 3.1415f;
+
+Glyph::Glyph(QQuickItem *parent)
+    : QQuickItem(parent)
+    , m_size(0)
+{
+    setFlag(QQuickItem::ItemHasContents);
+}
+
+void Glyph::setSize(qreal size)
+{
+    if (size == m_size)
+        return;
+
+    m_size = size;
+    setWidth(size);
+    setHeight(size);
+    emit sizeChanged();
+    update();
+}
+
+void Glyph::setColor(const QColor &color)
+{
+    if (color == m_color)
+        return;
+
+    m_color = color;
+    emit colorChanged();
+    update();
+}
+
+int calculateCircleVertexCount(qreal radius)
+{
+    // 10 * sqrt(r) \approx 2*pi / acos(1 - 1 / (4*r))
+    return (int) (10.0 * sqrt(radius));
+}
+
+void updateCircleGeometry(QSGGeometry *geometry, const QRectF &bounds)
+{
+    int vertexCount = geometry->vertexCount();
+    float cy = bounds.center().x();
+    float cx = bounds.center().y();
+
+    float theta = 2 * PI / float(vertexCount);
+    float c = cosf(theta);
+    float s = sinf(theta);
+    float x = bounds.width() / 2;
+    float y = 0;
+
+    QSGGeometry::Point2D *vertexData = geometry->vertexDataAsPoint2D();
+    for (int i = 0; i < vertexCount; i++) {
+        vertexData[i].set(x + cx, y + cy);
+
+        float t = x;
+        x = c*x - s*y;
+        y = s*t + c*y;
+    }
+}
+
+QSGNode *Glyph::updatePaintNode(QSGNode *oldNode, UpdatePaintNodeData *)
+{
+    QSGGeometryNode *node = 0;
+    QSGGeometry *geometry = 0;
+    QSGFlatColorMaterial *material = 0;
+    int vertexCount = calculateCircleVertexCount(m_size / 2);
+
+    if (!oldNode) {
+        node = new QSGGeometryNode;
+
+        geometry = new QSGGeometry(QSGGeometry::defaultAttributes_Point2D(), vertexCount);
+        geometry->setDrawingMode(GL_POLYGON);
+        node->setGeometry(geometry);
+        node->setFlag(QSGNode::OwnsGeometry);
+
+        material = new QSGFlatColorMaterial;
+        material->setColor(m_color);
+        node->setMaterial(material);
+        node->setFlag(QSGNode::OwnsMaterial);
+    } else {
+        node = static_cast<QSGGeometryNode *>(oldNode);
+        geometry = node->geometry();
+        geometry->allocate(vertexCount);
+    }
+
+    updateCircleGeometry(geometry, boundingRect());
+    node->markDirty(QSGNode::DirtyGeometry);
+
+    material = static_cast<QSGFlatColorMaterial *>(node->material());
+    if (material->color() != m_color) {
+        material->setColor(m_color);
+        node->markDirty(QSGNode::DirtyMaterial);
+    }
+
+    return node;
+}
diff --git a/glyph.h b/glyph.h
new file mode 100644
index 0000000..7532d04
--- /dev/null
+++ b/glyph.h
@@ -0,0 +1,42 @@
+#ifndef GLYPH_H
+#define GLYPH_H
+
+#include <QColor>
+#include <QtQuick/QQuickItem>
+
+class Glyph : public QQuickItem
+{
+    Q_OBJECT
+    Q_PROPERTY(QColor color READ color WRITE setColor NOTIFY colorChanged)
+    Q_PROPERTY(qreal size READ size WRITE setSize NOTIFY sizeChanged)
+
+public:
+    enum GlyphType {
+        GLYPH_CIRCLE,
+        GLYPH_SQUARE,
+        GLYPH_STAR,
+        GLYPH_CROSS
+    };
+
+    Glyph(QQuickItem *parent = 0);
+
+    QSGNode *updatePaintNode(QSGNode *oldNode, UpdatePaintNodeData *);
+
+    QColor color() const { return m_color; }
+    void setColor(const QColor &color);
+
+    qreal size() const { return m_size; }
+    void setSize(qreal size);
+
+signals:
+    void colorChanged();
+    void sizeChanged();
+
+public slots:
+
+private:
+    QColor m_color;
+    qreal m_size;
+};
+
+#endif // GLYPH_H
diff --git a/lamp.cpp b/lamp.cpp
new file mode 100644
index 0000000..a9d7a8b
--- /dev/null
+++ b/lamp.cpp
@@ -0,0 +1,50 @@
+#include "mp.h"
+
+#include <algorithm>
+
+arma::mat mp::lamp(const arma::mat &X, const arma::uvec &sampleIndices, const arma::mat &Ys)
+{
+    arma::mat Xs = X.rows(sampleIndices);
+    arma::uword sampleSize = sampleIndices.n_elem;
+    arma::mat projection(X.n_rows, 2);
+
+    for (arma::uword i = 0; i < X.n_rows; i++) {
+        arma::rowvec point = X.row(i);
+
+        // calculate alphas
+        arma::rowvec alphas(sampleSize);
+        for (arma::uword j = 0; j < sampleSize; j++) {
+            double dist = arma::accu(arma::square(Xs.row(j) - point));
+            alphas[j] = 1. / std::max(dist, mp::EPSILON);
+        }
+
+        double alphas_sum = arma::accu(alphas);
+        arma::rowvec alphas_sqrt = arma::sqrt(alphas);
+
+        // calculate \tilde{X} and \tilde{Y}
+        arma::rowvec Xtil = arma::sum(alphas * Xs, 0) / alphas_sum;
+        arma::rowvec Ytil = arma::sum(alphas * Ys, 0) / alphas_sum;
+
+        // calculate \hat{X} and \hat{Y}
+        arma::mat Xhat = Xs;
+        Xhat.each_row() -= Xtil;
+        arma::mat Yhat = Ys;
+        Yhat.each_row() -= Ytil;
+
+        // calculate A and B
+        arma::mat At = Xhat.t();
+        At.each_row() %= alphas_sqrt;
+        arma::mat B = Yhat;
+        B.each_col() %= alphas_sqrt.t();
+
+        arma::mat U, V;
+        arma::vec s;
+        arma::svd(U, s, V, At * B);
+        arma::mat M = U.cols(0, 1) * V.t();
+
+        // the projection of point i
+        projection.row(i) = (point - Xtil) * M + Ytil;
+    }
+
+    return projection;
+}
diff --git a/main.cpp b/main.cpp
new file mode 100644
index 0000000..863a20d
--- /dev/null
+++ b/main.cpp
@@ -0,0 +1,48 @@
+#include <cmath>
+#include <QGuiApplication>
+#include <QtQuick/QQuickView>
+
+#include "mp.h"
+#include "glyph.h"
+#include "scatterplot.h"
+
+arma::uvec getSample(arma::uword n)
+{
+    return arma::randi<arma::uvec>((arma::uword) 3*sqrt(n), arma::distr_param(0, n-1));
+}
+
+arma::mat getProjection(const arma::mat &X)
+{
+    arma::uword n = X.n_rows;
+    arma::uvec sampleIndices = getSample(n);
+    arma::mat Ys = arma::randn(sampleIndices.n_elem, 2);
+    Ys = mp::forceScheme(mp::dist(X.rows(sampleIndices)), Ys);
+    return mp::lamp(X, sampleIndices, Ys);
+}
+
+int main(int argc, char **argv)
+{
+    QGuiApplication app(argc, argv);
+
+    qmlRegisterType<Scatterplot>("PM", 1, 0, "Scatterplot");
+    qmlRegisterType<Glyph>("PM", 1, 0, "Glyph");
+
+    QQuickView view;
+    QSurfaceFormat format = view.format();
+    format.setSamples(16);
+    view.setFormat(format);
+    view.setResizeMode(QQuickView::SizeRootObjectToView);
+    view.setSource(QUrl("qrc:///main_view.qml"));
+
+    if (argc > 1) {
+        arma::mat X;
+        X.load(argv[1], arma::raw_ascii);
+        arma::mat projection = getProjection(X);
+        Scatterplot *plot = view.rootObject()->findChild<Scatterplot *>("plot", Qt::FindDirectChildrenOnly);
+        plot->setData(projection);
+    }
+
+    view.show();
+
+    return app.exec();
+}
diff --git a/main_view.qml b/main_view.qml
new file mode 100644
index 0000000..bf143fc
--- /dev/null
+++ b/main_view.qml
@@ -0,0 +1,40 @@
+import QtQuick 2.0
+import PM 1.0
+
+Item {
+    width: 480
+    height: 480
+
+    Scatterplot {
+        id: plot
+        objectName: "plot"
+        anchors.fill: parent
+    }
+
+    Item {
+        id: messageBox
+        anchors.right: parent.right
+        anchors.left: parent.left
+        anchors.bottom: parent.bottom
+
+        Rectangle {
+            color: Qt.rgba(1, 1, 1, 0.7)
+            radius: 5
+            border.width: 1
+            border.color: "white"
+            anchors.fill: messageLabel
+            anchors.margins: -10
+        }
+
+        Text {
+            id: messageLabel
+            color: "black"
+            wrapMode: Text.WordWrap
+            text: "The background here is a squircle rendered with raw OpenGL using the 'beforeRender()' signal in QQuickWindow. This text label and its border is rendered using QML"
+            anchors.right: parent.right
+            anchors.left: parent.left
+            anchors.bottom: parent.bottom
+            anchors.margins: 20
+        }
+    }
+}
diff --git a/mp.h b/mp.h
new file mode 100644
index 0000000..b3f9c75
--- /dev/null
+++ b/mp.h
@@ -0,0 +1,11 @@
+#include <armadillo>
+
+namespace mp {
+
+static const double EPSILON = 1e-3;
+double euclidean(const arma::rowvec &x1, const arma::rowvec &x2);
+arma::mat dist(const arma::mat &X, double (*distCalc)(const arma::rowvec &, const arma::rowvec &) = euclidean);
+arma::mat lamp(const arma::mat &X, const arma::uvec &sampleIndices, const arma::mat &Ys);
+arma::mat forceScheme(const arma::mat &D, arma::mat &Y, size_t maxIter = 20, double tol = 1e-3, double fraction = 0.25);
+
+} // namespace mp
diff --git a/pm.pro b/pm.pro
new file mode 100644
index 0000000..a8d87d7
--- /dev/null
+++ b/pm.pro
@@ -0,0 +1,16 @@
+QT += qml quick
+
+QMAKE_LIBS += -larmadillo
+HEADERS += glyph.h \
+    scatterplot.h \
+    mp.h
+SOURCES += main.cpp \
+    glyph.cpp \
+    scatterplot.cpp \
+    lamp.cpp \
+    forceScheme.cpp \
+    dist.cpp
+RESOURCES += pm.qrc
+
+target.path = .
+INSTALLS += target
diff --git a/pm.qrc b/pm.qrc
new file mode 100644
index 0000000..a3059be
--- /dev/null
+++ b/pm.qrc
@@ -0,0 +1,5 @@
+<RCC>
+    <qresource prefix="/">
+        <file>main_view.qml</file>
+    </qresource>
+</RCC>
diff --git a/scatterplot.cpp b/scatterplot.cpp
new file mode 100644
index 0000000..e775661
--- /dev/null
+++ b/scatterplot.cpp
@@ -0,0 +1,60 @@
+#include "scatterplot.h"
+
+#include <cstdio>
+#include "glyph.h"
+#include <QSGNode>
+#include <QSGGeometryNode>
+
+const int GLYPH_SIZE = 10;
+
+Scatterplot::Scatterplot()
+{
+    setFlag(QQuickItem::ItemHasContents);
+}
+
+Scatterplot::~Scatterplot()
+{
+}
+
+void Scatterplot::setData(const arma::mat &data)
+{
+    if (data.n_cols != 2)
+        return;
+
+    m_data = data;
+    qreal xmin = m_data.col(0).min(),
+          xmax = m_data.col(0).max(),
+          ymin = m_data.col(1).min(),
+          ymax = m_data.col(1).max();
+
+    for (arma::uword i = 0; i < m_data.n_rows; i++) {
+        arma::rowvec row = m_data.row(i);
+
+        Glyph *glyph = new Glyph();
+
+        glyph->setSize(5);
+        glyph->setX((row[0] - xmin) / (xmax - xmin) * width());
+        glyph->setY((row[1] - ymin) / (ymax - ymin) * height());
+
+        glyph->setParent(this);
+    }
+
+    update();
+}
+
+QSGNode *Scatterplot::updatePaintNode(QSGNode *oldNode, UpdatePaintNodeData *)
+{
+    QSGNode *node = 0;
+
+    if (!oldNode) {
+        node = new QSGNode;
+        for (QObjectList::const_iterator it = children().begin(); it != children().end(); it++)
+            node->appendChildNode(static_cast<Glyph *>(*it)->updatePaintNode(0, 0));
+    } else {
+        node = static_cast<QSGNode *>(oldNode);
+    }
+
+    node->markDirty(QSGNode::DirtyGeometry);
+
+    return node;
+}
diff --git a/scatterplot.h b/scatterplot.h
new file mode 100644
index 0000000..f38ce78
--- /dev/null
+++ b/scatterplot.h
@@ -0,0 +1,26 @@
+#ifndef SCATTERPLOT_H
+#define SCATTERPLOT_H
+
+#include <armadillo>
+#include <vector>
+#include <QQuickItem>
+
+class Scatterplot : public QQuickItem
+{
+    Q_OBJECT
+public:
+    Scatterplot();
+    ~Scatterplot();
+
+    QSGNode *updatePaintNode(QSGNode *oldNode, UpdatePaintNodeData *);
+    void setData(const arma::mat &data);
+
+signals:
+
+public slots:
+
+private:
+    arma::mat m_data;
+};
+
+#endif // SCATTERPLOT_H