QSimplex: Adding the Simplex solver for linear systems

This implementation is based on the iterative algorithm presented in
http://en.wikibooks.org/wiki/Operations_Research/The_Simplex_Method.
It is capable of giving us the solution to any n variable LP model.
Although we focused on QGraphicsAnchorLayout, the solver is general
enough for being reused in any other Qt classes if needed.

Signed-off-by: Anselmo Lacerda S. de Melo <anselmo.melo@openbossa.org>
Signed-off-by: Caio Marcelo de Oliveira Filho <caio.oliveira@openbossa.org>
Signed-off-by: Eduardo M. Fleury <eduardo.fleury@openbossa.org>
Signed-off-by: Jesus Sanchez-Palencia <jesus.palencia@openbossa.org>

3 files changed, 488 insertions, 2 deletions

diff --git a/src/gui/graphicsview/graphicsview.pri b/src/gui/graphicsview/graphicsview.pri
index ce35721..49ba70b 100644
--- a/src/gui/graphicsview/graphicsview.pri
+++ b/src/gui/graphicsview/graphicsview.pri
@@ -21,7 +21,8 @@ HEADERS += graphicsview/qgraphicsgridlayout.h \
     graphicsview/qgraphicswidget.h \
     graphicsview/qgraphicswidget_p.h \
     graphicsview/qgridlayoutengine_p.h \
-    graphicsview/qgraph_p.h
+    graphicsview/qgraph_p.h \
+    graphicsview/qsimplex_p.h
 
 SOURCES += graphicsview/qgraphicsgridlayout.cpp \
     graphicsview/qgraphicsitem.cpp \
@@ -40,4 +41,5 @@ SOURCES += graphicsview/qgraphicsgridlayout.cpp \
     graphicsview/qgraphicsview.cpp \
     graphicsview/qgraphicswidget.cpp \
     graphicsview/qgraphicswidget_p.cpp \
-    graphicsview/qgridlayoutengine.cpp
+    graphicsview/qgridlayoutengine.cpp \
+    graphicsview/qsimplex_p.cpp
diff --git a/src/gui/graphicsview/qsimplex_p.cpp b/src/gui/graphicsview/qsimplex_p.cpp
new file mode 100644
index 0000000..4ed11c4
--- /dev/null
+++ b/src/gui/graphicsview/qsimplex_p.cpp
@@ -0,0 +1,364 @@
+#include "qsimplex_p.h"
+
+#include <QtCore/qset.h>
+#include <QtCore/qdebug.h>
+
+#include <stdlib.h>
+
+QSimplex::QSimplex() : rows(0), columns(0), firstArtificial(0), matrix(0)
+{
+
+}
+
+QSimplex::~QSimplex()
+{
+    clearDataStructures();
+}
+
+void QSimplex::clearDataStructures()
+{
+    if (matrix == 0)
+        return;
+
+    // Matrix
+    rows = 0;
+    columns = 0;
+    firstArtificial = 0;
+    free(matrix);
+    matrix = 0;
+
+    // Constraints
+    for (int i = 0; i < constraints.size(); ++i) {
+        delete constraints[i]->helper.first;
+        constraints[i]->helper.first = 0;
+        constraints[i]->helper.second = 0.0;
+        delete constraints[i]->artificial;
+        constraints[i]->artificial = 0;
+    }
+    constraints.clear();
+
+    // Other
+    variables.clear();
+    objective = 0;
+}
+
+void QSimplex::setConstraints(const QList<QSimplexConstraint *> newConstraints)
+{
+    clearDataStructures();
+
+    if (newConstraints.isEmpty())
+        return;
+    constraints = newConstraints;
+
+    // Set Variables direct mapping
+    QSet<QSimplexVariable *> variablesSet;
+    for (int i = 0; i < constraints.size(); ++i)
+        variablesSet += \
+            QSet<QSimplexVariable *>::fromList(constraints[i]->variables.keys());
+    variables = variablesSet.toList();
+
+    // Set Variables reverse mapping
+    for (int i = 0; i < variables.size(); ++i) {
+        // The variable "0" goes at the column "1", etc...
+        variables[i]->index = i + 1;
+    }
+
+    // Normalize Constraints
+    int variableIndex = variables.size();
+    QList <QSimplexVariable *> artificialList;
+
+    for (int i = 0; i < constraints.size(); ++i) {
+        QSimplexVariable *slack;
+        QSimplexVariable *surplus;
+        QSimplexVariable *artificial;
+
+        Q_ASSERT(constraints[i]->helper.first == 0);
+        Q_ASSERT(constraints[i]->artificial == 0);
+
+        switch(constraints[i]->ratio) {
+        case QSimplexConstraint::LessOrEqual:
+            slack = new QSimplexVariable;
+            slack->index = ++variableIndex;
+            constraints[i]->helper.first = slack;
+            constraints[i]->helper.second = 1.0;
+            break;
+        case QSimplexConstraint::MoreOrEqual:
+            surplus = new QSimplexVariable;
+            surplus->index = ++variableIndex;
+            constraints[i]->helper.first = surplus;
+            constraints[i]->helper.second = -1.0;
+            // fall through
+        case QSimplexConstraint::Equal:
+            artificial = new QSimplexVariable;
+            constraints[i]->artificial = artificial;
+            artificialList += constraints[i]->artificial;
+            break;
+        }
+    }
+
+    firstArtificial = variableIndex + 1;
+    for (int i = 0; i < artificialList.size(); ++i)
+        artificialList[i]->index = ++variableIndex;
+    artificialList.clear();
+
+    // Matrix
+
+    // One for each variable plus the Basic and BFS columns (first and last)
+    columns = variableIndex + 2;
+    // One for each constraint plus the objective function
+    rows = constraints.size() + 1;
+
+    matrix = (qreal *)malloc(sizeof(qreal) * columns * rows);
+    if (!matrix) {
+        qWarning() << "QSimplex: Unable to allocate memory!";
+        return;
+    }
+    for (int i = columns * rows - 1; i >= 0; --i)
+        matrix[i] = 0.0;
+
+    // Fill Matrix
+    for (int i = 1; i <= constraints.size(); ++i) {
+        QSimplexConstraint *c = constraints[i - 1];
+
+        if (c->artificial) {
+            // Will use artificial basic variable
+            setValueAt(i, 0, c->artificial->index);
+            setValueAt(i, c->artificial->index, 1.0);
+
+            if (c->helper.second != 0.0) {
+                // Surplus variable
+                setValueAt(i, c->helper.first->index, c->helper.second);
+            }
+        } else {
+            // Slack is used as the basic variable
+            Q_ASSERT(c->helper.second == 1.0);
+            setValueAt(i, 0, c->helper.first->index);
+            setValueAt(i, c->helper.first->index, 1.0);
+        }
+
+        QHash<QSimplexVariable *, qreal>::const_iterator iter;
+        for (iter = c->variables.constBegin();
+             iter != c->variables.constEnd();
+             ++iter) {
+            setValueAt(i, iter.key()->index, iter.value());
+        }
+
+        setValueAt(i, columns - 1, c->constant);
+    }
+
+    // Set temporary objective: -1 * sum_of_artificial_vars
+    for (int j = firstArtificial; j < columns - 1; ++j)
+        setValueAt(0, j, 1.0);
+
+    // Maximize our objective (artificial vars go to zero)
+    solveMaxHelper();
+
+    if (valueAt(0, columns - 1) != 0.0) {
+        qWarning() << "QSimplex: No feaseable solution!";
+        clearDataStructures();
+        return;
+    }
+
+    // Remove artificial variables
+    clearColumns(firstArtificial, columns - 2);
+}
+
+void QSimplex::solveMaxHelper()
+{
+    reducedRowEchelon();
+    while (iterate()) ;
+}
+
+void QSimplex::setObjective(QSimplexConstraint *newObjective)
+{
+    objective = newObjective;
+}
+
+void QSimplex::clearRow(int rowIndex)
+{
+    qreal *item = matrix + rowIndex * columns;
+    for (int i = 0; i < columns; ++i)
+        item[i] = 0.0;
+}
+
+void QSimplex::clearColumns(int first, int last)
+{
+    for (int i = 0; i < rows; ++i) {
+        qreal *row = matrix + i * columns;
+        for (int j = first; j <= last; ++j)
+            row[j] = 0.0;
+    }
+}
+
+void QSimplex::dumpMatrix()
+{
+    printf("---- Simplex Matrix ----\n");
+
+    printf("       ");
+    for (int j = 0; j < columns; ++j)
+        printf("  <% 2d >", j);
+    printf("\n");
+
+    for (int i = 0; i < rows; ++i) {
+        printf("Row %2d:", i);
+
+        qreal *row = matrix + i * columns;
+        for (int j = 0; j < columns; ++j) {
+            printf("  % 2.2f", row[j]);
+        }
+        printf("\n");
+    }
+    printf("------------------------\n\n");
+}
+
+void QSimplex::combineRows(int toIndex, int fromIndex, qreal factor)
+{
+    if (!factor)
+        return;
+
+    qreal *from = matrix + fromIndex * columns;
+    qreal *to = matrix + toIndex * columns;
+
+    for (int j = 1; j < columns; ++j) {
+        to[j] += factor * from[j];
+
+        // ### Avoid Numerical errors
+        if (qAbs(to[j]) < 0.0000000001)
+            to[j] = 0.0;
+    }
+}
+
+int QSimplex::findPivotColumn()
+{
+    qreal min = 0;
+    int minIndex = -1;
+
+    for (int j = 0; j < columns-1; ++j) {
+        if (valueAt(0, j) < min) {
+            min = valueAt(0, j);
+            minIndex = j;
+        }
+    }
+
+    return minIndex;
+}
+
+int QSimplex::pivotRowForColumn(int column)
+{
+    qreal min = 999999999999.0; // ###
+    int minIndex = -1;
+
+    for (int i = 1; i < rows; ++i) {
+        qreal divisor = valueAt(i, column);
+        if (divisor <= 0)
+            continue;
+
+        qreal quotient = valueAt(i, columns - 1) / divisor;
+        if (quotient < min) {
+            min = quotient;
+            minIndex = i;
+        }
+    }
+
+    return minIndex;
+}
+
+void QSimplex::reducedRowEchelon()
+{
+    for (int i = 1; i < rows; ++i) {
+        int factorInObjectiveRow = valueAt(i, 0);
+        combineRows(0, i, -1 * valueAt(0, factorInObjectiveRow));
+    }
+}
+
+bool QSimplex::iterate()
+{
+    // Find Pivot column
+    int pivotColumn = findPivotColumn();
+    if (pivotColumn == -1)
+        return false;
+
+    // Find Pivot row for column
+    int pivotRow = pivotRowForColumn(pivotColumn);
+    if (pivotRow == -1) {
+        qWarning() << "QSimplex: Unbounded problem!";
+        return false;
+    }
+
+    // Normalize Pivot Row
+    qreal pivot = valueAt(pivotRow, pivotColumn);
+    if (pivot != 1.0)
+        combineRows(pivotRow, pivotRow, (1.0 - pivot) / pivot);
+
+    // Update other rows
+    for (int row=0; row < rows; ++row) {
+        if (row == pivotRow)
+            continue;
+
+        combineRows(row, pivotRow, -1 * valueAt(row, pivotColumn));
+    }
+
+    // Update first column
+    setValueAt(pivotRow, 0, pivotColumn);
+
+    //    dumpMatrix();
+    //    printf("------------ end of iteration --------------\n");
+    return true;
+}
+
+qreal QSimplex::solveMin()
+{
+    // Remove old objective
+    clearRow(0);
+
+    // Set new objective
+    QHash<QSimplexVariable *, qreal>::const_iterator iter;
+    for (iter = objective->variables.constBegin();
+         iter != objective->variables.constEnd();
+         ++iter) {
+        setValueAt(0, iter.key()->index, iter.value());
+    }
+
+    solveMaxHelper();
+    collectResults();
+
+    return -1 * valueAt(0, columns - 1);
+}
+
+qreal QSimplex::solveMax()
+{
+    // Remove old objective
+    clearRow(0);
+
+    // Set new objective
+    QHash<QSimplexVariable *, qreal>::const_iterator iter;
+    for (iter = objective->variables.constBegin();
+         iter != objective->variables.constEnd();
+         ++iter) {
+        setValueAt(0, iter.key()->index, -1 * iter.value());
+    }
+
+    solveMaxHelper();
+    collectResults();
+
+    return valueAt(0, columns - 1);
+}
+
+void QSimplex::collectResults()
+{
+    // All variables are zero unless overridden below.
+
+    // ### Is this really needed? Is there any chance that an
+    // important variable remains as non-basic at the end of simplex?
+    for (int i = 0; i < variables.size(); ++i)
+        variables[i]->result = 0;
+
+    // Basic variables
+    // Update the variable indicated in the first column with the value
+    // in the last column.
+    for (int i = 1; i < rows; ++i) {
+        int index = valueAt(i, 0) - 1;
+        if (index < variables.size())
+            variables[index]->result = valueAt(i, columns - 1);
+    }
+}
diff --git a/src/gui/graphicsview/qsimplex_p.h b/src/gui/graphicsview/qsimplex_p.h
new file mode 100644
index 0000000..3881893
--- /dev/null
+++ b/src/gui/graphicsview/qsimplex_p.h
@@ -0,0 +1,120 @@
+/****************************************************************************
+**
+** Copyright (C) 1992-$THISYEAR$ $TROLLTECH$. All rights reserved.
+**
+** This file is part of the $MODULE$ of the Qt Toolkit.
+**
+** $TROLLTECH_DUAL_LICENSE$
+**
+** This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+** WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+**
+****************************************************************************/
+
+#ifndef QSIMPLEX_P_H
+#define QSIMPLEX_P_H
+
+#include <QtCore/qhash.h>
+#include <QtCore/qpair.h>
+
+struct QSimplexVariable
+{
+    QSimplexVariable() : result(0), index(0) {};
+
+    qreal result;
+    uint index;
+};
+
+
+/*!
+  \internal
+
+  Representation of a LP constraint like:
+
+    (c1 * X1) + (c2 * X2) + ...  =  K
+                             or <=  K
+                             or >=  K
+
+    Where (ci, Xi) are the pairs in "variables" and K the real "constant".
+*/
+struct QSimplexConstraint
+{
+    QSimplexConstraint() : constant(0), ratio(Equal), artificial(0) {};
+
+    enum Ratio {
+        LessOrEqual = 0,
+        Equal,
+        MoreOrEqual
+    };
+
+    QHash<QSimplexVariable *, qreal> variables;
+    qreal constant;
+    Ratio ratio;
+
+    QPair<QSimplexVariable *, qreal> helper;
+    QSimplexVariable * artificial;
+};
+
+
+class QSimplex
+{
+public:
+    QSimplex();
+    virtual ~QSimplex();
+
+    qreal solveMin();
+    qreal solveMax();
+    QList<QSimplexVariable *> constraintsVariables();
+
+    void setConstraints(const QList<QSimplexConstraint *> constraints);
+    void setObjective(QSimplexConstraint *objective);
+
+    void dumpMatrix();
+
+private:
+    // Matrix handling
+    qreal valueAt(int row, int column);
+    void setValueAt(int row, int column, qreal value);
+    void clearRow(int rowIndex);
+    void clearColumns(int first, int last);
+    void combineRows(int toIndex, int fromIndex, qreal factor);
+
+    // Simplex
+    int findPivotColumn();
+    int pivotRowForColumn(int column);
+    void reducedRowEchelon();
+    bool iterate();
+
+    // Helpers
+    void clearDataStructures();
+    void solveMaxHelper();
+    void collectResults();
+
+    QList<QSimplexConstraint *> constraints;
+    QList<QSimplexVariable *> variables;
+    QSimplexConstraint *objective;
+
+    int rows;
+    int columns;
+    int firstArtificial;
+
+    qreal *matrix;
+};
+
+inline QList<QSimplexVariable *> QSimplex::constraintsVariables()
+{
+    return variables;
+}
+
+inline qreal QSimplex::valueAt(int rowIndex, int columnIndex)
+{
+    return matrix[rowIndex * columns + columnIndex];
+}
+
+inline void QSimplex::setValueAt(int rowIndex, int columnIndex, qreal value)
+{
+    matrix[rowIndex * columns + columnIndex] = value;
+}
+
+
+#endif