From 0b6df071d94ae8f7c9cdd3c96506a1420129e471 Mon Sep 17 00:00:00 2001 From: Samuel Fadel Date: Fri, 15 May 2015 18:10:52 -0300 Subject: Initial commit. ForceScheme seems bugged. --- dist.cpp | 20 +++++++++++ forceScheme.cpp | 45 +++++++++++++++++++++++++ glyph.cpp | 103 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ glyph.h | 42 +++++++++++++++++++++++ lamp.cpp | 50 +++++++++++++++++++++++++++ main.cpp | 48 ++++++++++++++++++++++++++ main_view.qml | 40 ++++++++++++++++++++++ mp.h | 11 ++++++ pm.pro | 16 +++++++++ pm.qrc | 5 +++ scatterplot.cpp | 60 +++++++++++++++++++++++++++++++++ scatterplot.h | 26 ++++++++++++++ 12 files changed, 466 insertions(+) create mode 100644 dist.cpp create mode 100644 forceScheme.cpp create mode 100644 glyph.cpp create mode 100644 glyph.h create mode 100644 lamp.cpp create mode 100644 main.cpp create mode 100644 main_view.qml create mode 100644 mp.h create mode 100644 pm.pro create mode 100644 pm.qrc create mode 100644 scatterplot.cpp create mode 100644 scatterplot.h 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 +#include + +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 +#include +#include +#include +#include + +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(oldNode); + geometry = node->geometry(); + geometry->allocate(vertexCount); + } + + updateCircleGeometry(geometry, boundingRect()); + node->markDirty(QSGNode::DirtyGeometry); + + material = static_cast(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 +#include + +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 + +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 +#include +#include + +#include "mp.h" +#include "glyph.h" +#include "scatterplot.h" + +arma::uvec getSample(arma::uword n) +{ + return arma::randi((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("PM", 1, 0, "Scatterplot"); + qmlRegisterType("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("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 + +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 @@ + + + main_view.qml + + 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 +#include "glyph.h" +#include +#include + +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(*it)->updatePaintNode(0, 0)); + } else { + node = static_cast(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 +#include +#include + +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 -- cgit v1.2.3