1 /****************************************************************************

     2 **

     3 ** Copyright (C) 2009 Nokia Corporation and/or its subsidiary(-ies).

     4 ** All rights reserved.

     5 ** Contact: Nokia Corporation (qt-info@nokia.com)

     6 **

     7 ** This file is part of the QtGui module of the Qt Toolkit.

     8 **

     9 ** $QT_BEGIN_LICENSE:LGPL$

    10 ** No Commercial Usage

    11 ** This file contains pre-release code and may not be distributed.

    12 ** You may use this file in accordance with the terms and conditions

    13 ** contained in the Technology Preview License Agreement accompanying

    14 ** this package.

    15 **

    16 ** GNU Lesser General Public License Usage

    17 ** Alternatively, this file may be used under the terms of the GNU Lesser

    18 ** General Public License version 2.1 as published by the Free Software

    19 ** Foundation and appearing in the file LICENSE.LGPL included in the

    20 ** packaging of this file.  Please review the following information to

    21 ** ensure the GNU Lesser General Public License version 2.1 requirements

    22 ** will be met: http://www.gnu.org/licenses/old-licenses/lgpl-2.1.html.

    23 **

    24 ** In addition, as a special exception, Nokia gives you certain additional

    25 ** rights.  These rights are described in the Nokia Qt LGPL Exception

    26 ** version 1.1, included in the file LGPL_EXCEPTION.txt in this package.

    27 **

    28 ** If you have questions regarding the use of this file, please contact

    29 ** Nokia at qt-info@nokia.com.

    30 **

    31 **

    32 **

    33 **

    34 **

    35 **

    36 **

    37 **

    38 ** $QT_END_LICENSE$

    39 **

    40 ****************************************************************************/

    41 

    42 #include "qsimplex_p.h"

    43 

    44 #include <QtCore/qset.h>

    45 #include <QtCore/qdebug.h>

    46 

    47 #include <stdlib.h>

    48 

    49 QT_BEGIN_NAMESPACE

    50 

    51 /*!

    52   \internal

    53   \class QSimplex

    54 

    55   The QSimplex class is a Linear Programming problem solver based on the two-phase

    56   simplex method.

    57 

    58   It takes a set of QSimplexConstraints as its restrictive constraints and an

    59   additional QSimplexConstraint as its objective function. Then methods to maximize

    60   and minimize the problem solution are provided.

    61 

    62   The two-phase simplex method is based on the following steps:

    63   First phase:

    64   1.a) Modify the original, complex, and possibly not feasible problem, into a new,

    65        easy to solve problem.

    66   1.b) Set as the objective of the new problem, a feasible solution for the original

    67        complex problem.

    68   1.c) Run simplex to optimize the modified problem and check whether a solution for

    69        the original problem exists.

    70 

    71   Second phase:

    72   2.a) Go back to the original problem with the feasibl (but not optimal) solution

    73        found in the first phase.

    74   2.b) Set the original objective.

    75   3.c) Run simplex to optimize the original problem towards its optimal solution.

    76 */

    77 

    78 /*!

    79   \internal

    80 */

    81 QSimplex::QSimplex() : objective(0), rows(0), columns(0), firstArtificial(0), matrix(0)

    82 {

    83 }

    84 

    85 /*!

    86   \internal

    87 */

    88 QSimplex::~QSimplex()

    89 {

    90     clearDataStructures();

    91 }

    92 

    93 /*!

    94   \internal

    95 */

    96 void QSimplex::clearDataStructures()

    97 {

    98     if (matrix == 0)

    99         return;

   100 

   101     // Matrix

   102     rows = 0;

   103     columns = 0;

   104     firstArtificial = 0;

   105     free(matrix);

   106     matrix = 0;

   107 

   108     // Constraints

   109     for (int i = 0; i < constraints.size(); ++i) {

   110         delete constraints[i]->helper.first;

   111         constraints[i]->helper.first = 0;

   112         constraints[i]->helper.second = 0.0;

   113         delete constraints[i]->artificial;

   114         constraints[i]->artificial = 0;

   115     }

   116     constraints.clear();

   117 

   118     // Other

   119     variables.clear();

   120     objective = 0;

   121 }

   122 

   123 /*!

   124   \internal

   125   Sets the new constraints in the simplex solver and returns whether the problem

   126   is feasible.

   127 

   128   This method sets the new constraints, normalizes them, creates the simplex matrix

   129   and runs the first simplex phase.

   130 */

   131 bool QSimplex::setConstraints(const QList<QSimplexConstraint *> newConstraints)

   132 {

   133     ////////////////////////////

   134     // Reset to initial state //

   135     ////////////////////////////

   136     clearDataStructures();

   137 

   138     if (newConstraints.isEmpty())

   139         return true;    // we are ok with no constraints

   140     constraints = newConstraints;

   141 

   142     ///////////////////////////////////////

   143     // Prepare variables and constraints //

   144     ///////////////////////////////////////

   145 

   146     // Set Variables direct mapping.

   147     // "variables" is a list that provides a stable, indexed list of all variables

   148     // used in this problem.

   149     QSet<QSimplexVariable *> variablesSet;

   150     for (int i = 0; i < constraints.size(); ++i)

   151         variablesSet += \

   152             QSet<QSimplexVariable *>::fromList(constraints[i]->variables.keys());

   153     variables = variablesSet.toList();

   154 

   155     // Set Variables reverse mapping

   156     // We also need to be able to find the index for a given variable, to do that

   157     // we store in each variable its index.

   158     for (int i = 0; i < variables.size(); ++i) {

   159         // The variable "0" goes at the column "1", etc...

   160         variables[i]->index = i + 1;

   161     }

   162 

   163     // Normalize Constraints

   164     // In this step, we prepare the constraints in two ways:

   165     // Firstly, we modify all constraints of type "LessOrEqual" or "MoreOrEqual"

   166     // by the adding slack or surplus variables and making them "Equal" constraints.

   167     // Secondly, we need every single constraint to have a direct, easy feasible

   168     // solution. Constraints that have slack variables are already easy to solve,

   169     // to all the others we add artificial variables.

   170     //

   171     // At the end we modify the constraints as follows:

   172     //  - LessOrEqual: SLACK variable is added.

   173     //  - Equal: ARTIFICIAL variable is added.

   174     //  - More or Equal: ARTIFICIAL and SURPLUS variables are added.

   175     int variableIndex = variables.size();

   176     QList <QSimplexVariable *> artificialList;

   177 

   178     for (int i = 0; i < constraints.size(); ++i) {

   179         QSimplexVariable *slack;

   180         QSimplexVariable *surplus;

   181         QSimplexVariable *artificial;

   182 

   183         Q_ASSERT(constraints[i]->helper.first == 0);

   184         Q_ASSERT(constraints[i]->artificial == 0);

   185 

   186         switch(constraints[i]->ratio) {

   187         case QSimplexConstraint::LessOrEqual:

   188             slack = new QSimplexVariable;

   189             slack->index = ++variableIndex;

   190             constraints[i]->helper.first = slack;

   191             constraints[i]->helper.second = 1.0;

   192             break;

   193         case QSimplexConstraint::MoreOrEqual:

   194             surplus = new QSimplexVariable;

   195             surplus->index = ++variableIndex;

   196             constraints[i]->helper.first = surplus;

   197             constraints[i]->helper.second = -1.0;

   198             // fall through

   199         case QSimplexConstraint::Equal:

   200             artificial = new QSimplexVariable;

   201             constraints[i]->artificial = artificial;

   202             artificialList += constraints[i]->artificial;

   203             break;

   204         }

   205     }

   206 

   207     // All original, slack and surplus have already had its index set

   208     // at this point. We now set the index of the artificial variables

   209     // as to ensure they are at the end of the variable list and therefore

   210     // can be easily removed at the end of this method.

   211     firstArtificial = variableIndex + 1;

   212     for (int i = 0; i < artificialList.size(); ++i)

   213         artificialList[i]->index = ++variableIndex;

   214     artificialList.clear();

   215 

   216     /////////////////////////////

   217     // Fill the Simplex matrix //

   218     /////////////////////////////

   219 

   220     // One for each variable plus the Basic and BFS columns (first and last)

   221     columns = variableIndex + 2;

   222     // One for each constraint plus the objective function

   223     rows = constraints.size() + 1;

   224 

   225     matrix = (qreal *)malloc(sizeof(qreal) * columns * rows);

   226     if (!matrix) {

   227         qWarning() << "QSimplex: Unable to allocate memory!";

   228         return false;

   229     }

   230     for (int i = columns * rows - 1; i >= 0; --i)

   231         matrix[i] = 0.0;

   232 

   233     // Fill Matrix

   234     for (int i = 1; i <= constraints.size(); ++i) {

   235         QSimplexConstraint *c = constraints[i - 1];

   236 

   237         if (c->artificial) {

   238             // Will use artificial basic variable

   239             setValueAt(i, 0, c->artificial->index);

   240             setValueAt(i, c->artificial->index, 1.0);

   241 

   242             if (c->helper.second != 0.0) {

   243                 // Surplus variable

   244                 setValueAt(i, c->helper.first->index, c->helper.second);

   245             }

   246         } else {

   247             // Slack is used as the basic variable

   248             Q_ASSERT(c->helper.second == 1.0);

   249             setValueAt(i, 0, c->helper.first->index);

   250             setValueAt(i, c->helper.first->index, 1.0);

   251         }

   252 

   253         QHash<QSimplexVariable *, qreal>::const_iterator iter;

   254         for (iter = c->variables.constBegin();

   255              iter != c->variables.constEnd();

   256              ++iter) {

   257             setValueAt(i, iter.key()->index, iter.value());

   258         }

   259 

   260         setValueAt(i, columns - 1, c->constant);

   261     }

   262 

   263     // Set objective for the first-phase Simplex.

   264     // Z = -1 * sum_of_artificial_vars

   265     for (int j = firstArtificial; j < columns - 1; ++j)

   266         setValueAt(0, j, 1.0);

   267 

   268     // Maximize our objective (artificial vars go to zero)

   269     solveMaxHelper();

   270 

   271     // If there is a solution where the sum of all artificial

   272     // variables is zero, then all of them can be removed and yet

   273     // we will have a feasible (but not optimal) solution for the

   274     // original problem.

   275     // Otherwise, we clean up our structures and report there is

   276     // no feasible solution.

   277     if (valueAt(0, columns - 1) != 0.0) {

   278         qWarning() << "QSimplex: No feasible solution!";

   279         clearDataStructures();

   280         return false;

   281     }

   282 

   283     // Remove artificial variables. We already have a feasible

   284     // solution for the first problem, thus we don't need them

   285     // anymore.

   286     clearColumns(firstArtificial, columns - 2);

   287 

   288     return true;

   289 }

   290 

   291 /*!

   292   \internal

   293 

   294   Run simplex on the current matrix with the current objective.

   295 

   296   This is the iterative method. The matrix lines are combined

   297   as to modify the variable values towards the best solution possible.

   298   The method returns when the matrix is in the optimal state.

   299 */

   300 void QSimplex::solveMaxHelper()

   301 {

   302     reducedRowEchelon();

   303     while (iterate()) ;

   304 }

   305 

   306 /*!

   307   \internal

   308 */

   309 void QSimplex::setObjective(QSimplexConstraint *newObjective)

   310 {

   311     objective = newObjective;

   312 }

   313 

   314 /*!

   315   \internal

   316 */

   317 void QSimplex::clearRow(int rowIndex)

   318 {

   319     qreal *item = matrix + rowIndex * columns;

   320     for (int i = 0; i < columns; ++i)

   321         item[i] = 0.0;

   322 }

   323 

   324 /*!

   325   \internal

   326 */

   327 void QSimplex::clearColumns(int first, int last)

   328 {

   329     for (int i = 0; i < rows; ++i) {

   330         qreal *row = matrix + i * columns;

   331         for (int j = first; j <= last; ++j)

   332             row[j] = 0.0;

   333     }

   334 }

   335 

   336 /*!

   337   \internal

   338 */

   339 void QSimplex::dumpMatrix()

   340 {

   341     qDebug("---- Simplex Matrix ----\n");

   342 

   343     QString str(QLatin1String("       "));

   344     for (int j = 0; j < columns; ++j)

   345         str += QString::fromAscii("  <%1 >").arg(j, 2);

   346     qDebug("%s", qPrintable(str));

   347     for (int i = 0; i < rows; ++i) {

   348         str = QString::fromAscii("Row %1:").arg(i, 2);

   349 

   350         qreal *row = matrix + i * columns;

   351         for (int j = 0; j < columns; ++j)

   352             str += QString::fromAscii("%1").arg(row[j], 7, 'f', 2);

   353         qDebug("%s", qPrintable(str));

   354     }

   355     qDebug("------------------------\n");

   356 }

   357 

   358 /*!

   359   \internal

   360 */

   361 void QSimplex::combineRows(int toIndex, int fromIndex, qreal factor)

   362 {

   363     if (!factor)

   364         return;

   365 

   366     qreal *from = matrix + fromIndex * columns;

   367     qreal *to = matrix + toIndex * columns;

   368 

   369     for (int j = 1; j < columns; ++j) {

   370         qreal value = from[j];

   371 

   372         // skip to[j] = to[j] + factor*0.0

   373         if (value == 0.0)

   374             continue;

   375 

   376         to[j] += factor * value;

   377 

   378         // ### Avoid Numerical errors

   379         if (qAbs(to[j]) < 0.0000000001)

   380             to[j] = 0.0;

   381     }

   382 }

   383 

   384 /*!

   385   \internal

   386 */

   387 int QSimplex::findPivotColumn()

   388 {

   389     qreal min = 0;

   390     int minIndex = -1;

   391 

   392     for (int j = 0; j < columns-1; ++j) {

   393         if (valueAt(0, j) < min) {

   394             min = valueAt(0, j);

   395             minIndex = j;

   396         }

   397     }

   398 

   399     return minIndex;

   400 }

   401 

   402 /*!

   403   \internal

   404 

   405   For a given pivot column, find the pivot row. That is, the row with the

   406   minimum associated "quotient" where:

   407 

   408   - quotient is the division of the value in the last column by the value

   409     in the pivot column.

   410   - rows with value less or equal to zero are ignored

   411   - if two rows have the same quotient, lines are chosen based on the

   412     highest variable index (value in the first column)

   413 

   414   The last condition avoids a bug where artificial variables would be

   415   left behind for the second-phase simplex, and with 'good'

   416   constraints would be removed before it, what would lead to incorrect

   417   results.

   418 */

   419 int QSimplex::pivotRowForColumn(int column)

   420 {

   421     qreal min = qreal(999999999999.0); // ###

   422     int minIndex = -1;

   423 

   424     for (int i = 1; i < rows; ++i) {

   425         qreal divisor = valueAt(i, column);

   426         if (divisor <= 0)

   427             continue;

   428 

   429         qreal quotient = valueAt(i, columns - 1) / divisor;

   430         if (quotient < min) {

   431             min = quotient;

   432             minIndex = i;

   433         } else if ((quotient == min) && (valueAt(i, 0) > valueAt(minIndex, 0))) {

   434             minIndex = i;

   435         }

   436     }

   437 

   438     return minIndex;

   439 }

   440 

   441 /*!

   442   \internal

   443 */

   444 void QSimplex::reducedRowEchelon()

   445 {

   446     for (int i = 1; i < rows; ++i) {

   447         int factorInObjectiveRow = valueAt(i, 0);

   448         combineRows(0, i, -1 * valueAt(0, factorInObjectiveRow));

   449     }

   450 }

   451 

   452 /*!

   453   \internal

   454 

   455   Does one iteration towards a better solution for the problem.

   456   See 'solveMaxHelper'.

   457 */

   458 bool QSimplex::iterate()

   459 {

   460     // Find Pivot column

   461     int pivotColumn = findPivotColumn();

   462     if (pivotColumn == -1)

   463         return false;

   464 

   465     // Find Pivot row for column

   466     int pivotRow = pivotRowForColumn(pivotColumn);

   467     if (pivotRow == -1) {

   468         qWarning() << "QSimplex: Unbounded problem!";

   469         return false;

   470     }

   471 

   472     // Normalize Pivot Row

   473     qreal pivot = valueAt(pivotRow, pivotColumn);

   474     if (pivot != 1.0)

   475         combineRows(pivotRow, pivotRow, (1.0 - pivot) / pivot);

   476 

   477     // Update other rows

   478     for (int row=0; row < rows; ++row) {

   479         if (row == pivotRow)

   480             continue;

   481 

   482         combineRows(row, pivotRow, -1 * valueAt(row, pivotColumn));

   483     }

   484 

   485     // Update first column

   486     setValueAt(pivotRow, 0, pivotColumn);

   487 

   488     //    dumpMatrix();

   489     //    qDebug("------------ end of iteration --------------\n");

   490     return true;

   491 }

   492 

   493 /*!

   494   \internal

   495 

   496   Both solveMin and solveMax are interfaces to this method.

   497 

   498   The enum solverFactor admits 2 values: Minimum (-1) and Maximum (+1).

   499 

   500   This method sets the original objective and runs the second phase

   501   Simplex to obtain the optimal solution for the problem. As the internal

   502   simplex solver is only able to _maximize_ objectives, we handle the

   503   minimization case by inverting the original objective and then

   504   maximizing it.

   505 */

   506 qreal QSimplex::solver(solverFactor factor)

   507 {

   508     // Remove old objective

   509     clearRow(0);

   510 

   511     // Set new objective

   512     QHash<QSimplexVariable *, qreal>::const_iterator iter;

   513     for (iter = objective->variables.constBegin();

   514          iter != objective->variables.constEnd();

   515          ++iter) {

   516         setValueAt(0, iter.key()->index, -1 * factor * iter.value());

   517     }

   518 

   519     solveMaxHelper();

   520     collectResults();

   521 

   522 #ifdef QT_DEBUG

   523     for (int i = 0; i < constraints.size(); ++i) {

   524         Q_ASSERT(constraints[i]->isSatisfied());

   525     }

   526 #endif

   527 

   528     return factor * valueAt(0, columns - 1);

   529 }

   530 

   531 /*!

   532   \internal

   533   Minimize the original objective.

   534 */

   535 qreal QSimplex::solveMin()

   536 {

   537     return solver(Minimum);

   538 }

   539 

   540 /*!

   541   \internal

   542   Maximize the original objective.

   543 */

   544 qreal QSimplex::solveMax()

   545 {

   546     return solver(Maximum);

   547 }

   548 

   549 /*!

   550   \internal

   551 

   552   Reads results from the simplified matrix and saves them in the

   553   "result" member of each QSimplexVariable.

   554 */

   555 void QSimplex::collectResults()

   556 {

   557     // All variables are zero unless overridden below.

   558 

   559     // ### Is this really needed? Is there any chance that an

   560     // important variable remains as non-basic at the end of simplex?

   561     for (int i = 0; i < variables.size(); ++i)

   562         variables[i]->result = 0;

   563 

   564     // Basic variables

   565     // Update the variable indicated in the first column with the value

   566     // in the last column.

   567     for (int i = 1; i < rows; ++i) {

   568         int index = valueAt(i, 0) - 1;

   569         if (index < variables.size())

   570             variables[index]->result = valueAt(i, columns - 1);

   571     }

   572 }

   573 

   574 QT_END_NAMESPACE