1/****************************************************************************
2**
3** Copyright (C) 2016 The Qt Company Ltd.
4** Contact: https://www.qt.io/licensing/
5**
6** This file is part of the QtWidgets module of the Qt Toolkit.
7**
8** $QT_BEGIN_LICENSE:LGPL$
9** Commercial License Usage
10** Licensees holding valid commercial Qt licenses may use this file in
11** accordance with the commercial license agreement provided with the
12** Software or, alternatively, in accordance with the terms contained in
13** a written agreement between you and The Qt Company. For licensing terms
14** and conditions see https://www.qt.io/terms-conditions. For further
15** information use the contact form at https://www.qt.io/contact-us.
16**
17** GNU Lesser General Public License Usage
18** Alternatively, this file may be used under the terms of the GNU Lesser
19** General Public License version 3 as published by the Free Software
20** Foundation and appearing in the file LICENSE.LGPL3 included in the
21** packaging of this file. Please review the following information to
22** ensure the GNU Lesser General Public License version 3 requirements
23** will be met: https://www.gnu.org/licenses/lgpl-3.0.html.
24**
25** GNU General Public License Usage
26** Alternatively, this file may be used under the terms of the GNU
27** General Public License version 2.0 or (at your option) the GNU General
28** Public license version 3 or any later version approved by the KDE Free
29** Qt Foundation. The licenses are as published by the Free Software
30** Foundation and appearing in the file LICENSE.GPL2 and LICENSE.GPL3
31** included in the packaging of this file. Please review the following
32** information to ensure the GNU General Public License requirements will
33** be met: https://www.gnu.org/licenses/gpl-2.0.html and
34** https://www.gnu.org/licenses/gpl-3.0.html.
35**
36** $QT_END_LICENSE$
37**
38****************************************************************************/
39
40#include "qsimplex_p.h"
41
42#include <QtCore/qset.h>
43#include <QtCore/qdebug.h>
44
45#include <stdlib.h>
46
47QT_BEGIN_NAMESPACE
48
49/*!
50 \internal
51 \class QSimplex
52
53 The QSimplex class is a Linear Programming problem solver based on the two-phase
54 simplex method.
55
56 It takes a set of QSimplexConstraints as its restrictive constraints and an
57 additional QSimplexConstraint as its objective function. Then methods to maximize
58 and minimize the problem solution are provided.
59
60 The two-phase simplex method is based on the following steps:
61 First phase:
62 1.a) Modify the original, complex, and possibly not feasible problem, into a new,
63 easy to solve problem.
64 1.b) Set as the objective of the new problem, a feasible solution for the original
65 complex problem.
66 1.c) Run simplex to optimize the modified problem and check whether a solution for
67 the original problem exists.
68
69 Second phase:
70 2.a) Go back to the original problem with the feasibl (but not optimal) solution
71 found in the first phase.
72 2.b) Set the original objective.
73 3.c) Run simplex to optimize the original problem towards its optimal solution.
74*/
75
76/*!
77 \internal
78*/
79QSimplex::QSimplex() : objective(nullptr), rows(0), columns(0), firstArtificial(0), matrix(nullptr)
80{
81}
82
83/*!
84 \internal
85*/
86QSimplex::~QSimplex()
87{
88 clearDataStructures();
89}
90
91/*!
92 \internal
93*/
94void QSimplex::clearDataStructures()
95{
96 if (matrix == nullptr)
97 return;
98
99 // Matrix
100 rows = 0;
101 columns = 0;
102 firstArtificial = 0;
103 free(matrix);
104 matrix = nullptr;
105
106 // Constraints
107 for (int i = 0; i < constraints.size(); ++i) {
108 delete constraints[i]->helper.first;
109 delete constraints[i]->artificial;
110 delete constraints[i];
111 }
112 constraints.clear();
113
114 // Other
115 variables.clear();
116 objective = nullptr;
117}
118
119/*!
120 \internal
121 Sets the new constraints in the simplex solver and returns whether the problem
122 is feasible.
123
124 This method sets the new constraints, normalizes them, creates the simplex matrix
125 and runs the first simplex phase.
126*/
127bool QSimplex::setConstraints(const QList<QSimplexConstraint *> &newConstraints)
128{
129 ////////////////////////////
130 // Reset to initial state //
131 ////////////////////////////
132 clearDataStructures();
133
134 if (newConstraints.isEmpty())
135 return true; // we are ok with no constraints
136
137 // Make deep copy of constraints. We need this copy because we may change
138 // them in the simplification method.
139 for (int i = 0; i < newConstraints.size(); ++i) {
140 QSimplexConstraint *c = new QSimplexConstraint;
141 c->constant = newConstraints[i]->constant;
142 c->ratio = newConstraints[i]->ratio;
143 c->variables = newConstraints[i]->variables;
144 constraints << c;
145 }
146
147 // Remove constraints of type Var == K and replace them for their value.
148 if (!simplifyConstraints(&constraints)) {
149 qWarning("QSimplex: No feasible solution!");
150 clearDataStructures();
151 return false;
152 }
153
154 ///////////////////////////////////////
155 // Prepare variables and constraints //
156 ///////////////////////////////////////
157
158 // Set Variables direct mapping.
159 // "variables" is a list that provides a stable, indexed list of all variables
160 // used in this problem.
161 QSet<QSimplexVariable *> variablesSet;
162 for (int i = 0; i < constraints.size(); ++i) {
163 const auto &v = constraints.at(i)->variables;
164 for (auto it = v.cbegin(), end = v.cend(); it != end; ++it)
165 variablesSet.insert(it.key());
166 }
167 variables = variablesSet.values();
168
169 // Set Variables reverse mapping
170 // We also need to be able to find the index for a given variable, to do that
171 // we store in each variable its index.
172 for (int i = 0; i < variables.size(); ++i) {
173 // The variable "0" goes at the column "1", etc...
174 variables[i]->index = i + 1;
175 }
176
177 // Normalize Constraints
178 // In this step, we prepare the constraints in two ways:
179 // Firstly, we modify all constraints of type "LessOrEqual" or "MoreOrEqual"
180 // by the adding slack or surplus variables and making them "Equal" constraints.
181 // Secondly, we need every single constraint to have a direct, easy feasible
182 // solution. Constraints that have slack variables are already easy to solve,
183 // to all the others we add artificial variables.
184 //
185 // At the end we modify the constraints as follows:
186 // - LessOrEqual: SLACK variable is added.
187 // - Equal: ARTIFICIAL variable is added.
188 // - More or Equal: ARTIFICIAL and SURPLUS variables are added.
189 int variableIndex = variables.size();
190 QList <QSimplexVariable *> artificialList;
191
192 for (int i = 0; i < constraints.size(); ++i) {
193 QSimplexVariable *slack;
194 QSimplexVariable *surplus;
195 QSimplexVariable *artificial;
196
197 Q_ASSERT(constraints[i]->helper.first == 0);
198 Q_ASSERT(constraints[i]->artificial == nullptr);
199
200 switch(constraints[i]->ratio) {
201 case QSimplexConstraint::LessOrEqual:
202 slack = new QSimplexVariable;
203 slack->index = ++variableIndex;
204 constraints[i]->helper.first = slack;
205 constraints[i]->helper.second = 1.0;
206 break;
207 case QSimplexConstraint::MoreOrEqual:
208 surplus = new QSimplexVariable;
209 surplus->index = ++variableIndex;
210 constraints[i]->helper.first = surplus;
211 constraints[i]->helper.second = -1.0;
212 Q_FALLTHROUGH();
213 case QSimplexConstraint::Equal:
214 artificial = new QSimplexVariable;
215 constraints[i]->artificial = artificial;
216 artificialList += constraints[i]->artificial;
217 break;
218 }
219 }
220
221 // All original, slack and surplus have already had its index set
222 // at this point. We now set the index of the artificial variables
223 // as to ensure they are at the end of the variable list and therefore
224 // can be easily removed at the end of this method.
225 firstArtificial = variableIndex + 1;
226 for (int i = 0; i < artificialList.size(); ++i)
227 artificialList[i]->index = ++variableIndex;
228 artificialList.clear();
229
230 /////////////////////////////
231 // Fill the Simplex matrix //
232 /////////////////////////////
233
234 // One for each variable plus the Basic and BFS columns (first and last)
235 columns = variableIndex + 2;
236 // One for each constraint plus the objective function
237 rows = constraints.size() + 1;
238
239 matrix = (qreal *)malloc(sizeof(qreal) * columns * rows);
240 if (!matrix) {
241 qWarning("QSimplex: Unable to allocate memory!");
242 return false;
243 }
244 for (int i = columns * rows - 1; i >= 0; --i)
245 matrix[i] = 0.0;
246
247 // Fill Matrix
248 for (int i = 1; i <= constraints.size(); ++i) {
249 QSimplexConstraint *c = constraints[i - 1];
250
251 if (c->artificial) {
252 // Will use artificial basic variable
253 setValueAt(i, 0, c->artificial->index);
254 setValueAt(i, c->artificial->index, 1.0);
255
256 if (c->helper.second != 0.0) {
257 // Surplus variable
258 setValueAt(i, c->helper.first->index, c->helper.second);
259 }
260 } else {
261 // Slack is used as the basic variable
262 Q_ASSERT(c->helper.second == 1.0);
263 setValueAt(i, 0, c->helper.first->index);
264 setValueAt(i, c->helper.first->index, 1.0);
265 }
266
267 QHash<QSimplexVariable *, qreal>::const_iterator iter;
268 for (iter = c->variables.constBegin();
269 iter != c->variables.constEnd();
270 ++iter) {
271 setValueAt(i, iter.key()->index, iter.value());
272 }
273
274 setValueAt(i, columns - 1, c->constant);
275 }
276
277 // Set objective for the first-phase Simplex.
278 // Z = -1 * sum_of_artificial_vars
279 for (int j = firstArtificial; j < columns - 1; ++j)
280 setValueAt(0, j, 1.0);
281
282 // Maximize our objective (artificial vars go to zero)
283 solveMaxHelper();
284
285 // If there is a solution where the sum of all artificial
286 // variables is zero, then all of them can be removed and yet
287 // we will have a feasible (but not optimal) solution for the
288 // original problem.
289 // Otherwise, we clean up our structures and report there is
290 // no feasible solution.
291 if ((valueAt(0, columns - 1) != 0.0) && (qAbs(valueAt(0, columns - 1)) > 0.00001)) {
292 qWarning("QSimplex: No feasible solution!");
293 clearDataStructures();
294 return false;
295 }
296
297 // Remove artificial variables. We already have a feasible
298 // solution for the first problem, thus we don't need them
299 // anymore.
300 clearColumns(firstArtificial, columns - 2);
301
302 return true;
303}
304
305/*!
306 \internal
307
308 Run simplex on the current matrix with the current objective.
309
310 This is the iterative method. The matrix lines are combined
311 as to modify the variable values towards the best solution possible.
312 The method returns when the matrix is in the optimal state.
313*/
314void QSimplex::solveMaxHelper()
315{
316 reducedRowEchelon();
317 while (iterate()) ;
318}
319
320/*!
321 \internal
322*/
323void QSimplex::setObjective(QSimplexConstraint *newObjective)
324{
325 objective = newObjective;
326}
327
328/*!
329 \internal
330*/
331void QSimplex::clearRow(int rowIndex)
332{
333 qreal *item = matrix + rowIndex * columns;
334 for (int i = 0; i < columns; ++i)
335 item[i] = 0.0;
336}
337
338/*!
339 \internal
340*/
341void QSimplex::clearColumns(int first, int last)
342{
343 for (int i = 0; i < rows; ++i) {
344 qreal *row = matrix + i * columns;
345 for (int j = first; j <= last; ++j)
346 row[j] = 0.0;
347 }
348}
349
350/*!
351 \internal
352*/
353void QSimplex::dumpMatrix()
354{
355 qDebug("---- Simplex Matrix ----\n");
356
357 QString str(QLatin1String(" "));
358 for (int j = 0; j < columns; ++j)
359 str += QString::fromLatin1(" <%1 >").arg(j, 2);
360 qDebug("%s", qPrintable(str));
361 for (int i = 0; i < rows; ++i) {
362 str = QString::fromLatin1("Row %1:").arg(i, 2);
363
364 qreal *row = matrix + i * columns;
365 for (int j = 0; j < columns; ++j)
366 str += QString::fromLatin1("%1").arg(row[j], 7, 'f', 2);
367 qDebug("%s", qPrintable(str));
368 }
369 qDebug("------------------------\n");
370}
371
372/*!
373 \internal
374*/
375void QSimplex::combineRows(int toIndex, int fromIndex, qreal factor)
376{
377 if (!factor)
378 return;
379
380 qreal *from = matrix + fromIndex * columns;
381 qreal *to = matrix + toIndex * columns;
382
383 for (int j = 1; j < columns; ++j) {
384 qreal value = from[j];
385
386 // skip to[j] = to[j] + factor*0.0
387 if (value == 0.0)
388 continue;
389
390 to[j] += factor * value;
391
392 // ### Avoid Numerical errors
393 if (qAbs(to[j]) < 0.0000000001)
394 to[j] = 0.0;
395 }
396}
397
398/*!
399 \internal
400*/
401int QSimplex::findPivotColumn()
402{
403 qreal min = 0;
404 int minIndex = -1;
405
406 for (int j = 0; j < columns-1; ++j) {
407 if (valueAt(0, j) < min) {
408 min = valueAt(0, j);
409 minIndex = j;
410 }
411 }
412
413 return minIndex;
414}
415
416/*!
417 \internal
418
419 For a given pivot column, find the pivot row. That is, the row with the
420 minimum associated "quotient" where:
421
422 - quotient is the division of the value in the last column by the value
423 in the pivot column.
424 - rows with value less or equal to zero are ignored
425 - if two rows have the same quotient, lines are chosen based on the
426 highest variable index (value in the first column)
427
428 The last condition avoids a bug where artificial variables would be
429 left behind for the second-phase simplex, and with 'good'
430 constraints would be removed before it, what would lead to incorrect
431 results.
432*/
433int QSimplex::pivotRowForColumn(int column)
434{
435 qreal min = qreal(999999999999.0); // ###
436 int minIndex = -1;
437
438 for (int i = 1; i < rows; ++i) {
439 qreal divisor = valueAt(i, column);
440 if (divisor <= 0)
441 continue;
442
443 qreal quotient = valueAt(i, columns - 1) / divisor;
444 if (quotient < min) {
445 min = quotient;
446 minIndex = i;
447 } else if ((quotient == min) && (valueAt(i, 0) > valueAt(minIndex, 0))) {
448 minIndex = i;
449 }
450 }
451
452 return minIndex;
453}
454
455/*!
456 \internal
457*/
458void QSimplex::reducedRowEchelon()
459{
460 for (int i = 1; i < rows; ++i) {
461 int factorInObjectiveRow = valueAt(i, 0);
462 combineRows(0, i, -1 * valueAt(0, factorInObjectiveRow));
463 }
464}
465
466/*!
467 \internal
468
469 Does one iteration towards a better solution for the problem.
470 See 'solveMaxHelper'.
471*/
472bool QSimplex::iterate()
473{
474 // Find Pivot column
475 int pivotColumn = findPivotColumn();
476 if (pivotColumn == -1)
477 return false;
478
479 // Find Pivot row for column
480 int pivotRow = pivotRowForColumn(pivotColumn);
481 if (pivotRow == -1) {
482 qWarning("QSimplex: Unbounded problem!");
483 return false;
484 }
485
486 // Normalize Pivot Row
487 qreal pivot = valueAt(pivotRow, pivotColumn);
488 if (pivot != 1.0)
489 combineRows(pivotRow, pivotRow, (1.0 - pivot) / pivot);
490
491 // Update other rows
492 for (int row=0; row < rows; ++row) {
493 if (row == pivotRow)
494 continue;
495
496 combineRows(row, pivotRow, -1 * valueAt(row, pivotColumn));
497 }
498
499 // Update first column
500 setValueAt(pivotRow, 0, pivotColumn);
501
502 // dumpMatrix();
503 // qDebug("------------ end of iteration --------------\n");
504 return true;
505}
506
507/*!
508 \internal
509
510 Both solveMin and solveMax are interfaces to this method.
511
512 The enum SolverFactor admits 2 values: Minimum (-1) and Maximum (+1).
513
514 This method sets the original objective and runs the second phase
515 Simplex to obtain the optimal solution for the problem. As the internal
516 simplex solver is only able to _maximize_ objectives, we handle the
517 minimization case by inverting the original objective and then
518 maximizing it.
519*/
520qreal QSimplex::solver(SolverFactor factor)
521{
522 // Remove old objective
523 clearRow(0);
524
525 // Set new objective in the first row of the simplex matrix
526 qreal resultOffset = 0;
527 QHash<QSimplexVariable *, qreal>::const_iterator iter;
528 for (iter = objective->variables.constBegin();
529 iter != objective->variables.constEnd();
530 ++iter) {
531
532 // Check if the variable was removed in the simplification process.
533 // If so, we save its offset to the objective function and skip adding
534 // it to the matrix.
535 if (iter.key()->index == -1) {
536 resultOffset += iter.value() * iter.key()->result;
537 continue;
538 }
539
540 setValueAt(0, iter.key()->index, -1 * factor * iter.value());
541 }
542
543 solveMaxHelper();
544 collectResults();
545
546#ifdef QT_DEBUG
547 for (int i = 0; i < constraints.size(); ++i) {
548 Q_ASSERT(constraints[i]->isSatisfied());
549 }
550#endif
551
552 // Return the value calculated by the simplex plus the value of the
553 // fixed variables.
554 return (factor * valueAt(0, columns - 1)) + resultOffset;
555}
556
557/*!
558 \internal
559 Minimize the original objective.
560*/
561qreal QSimplex::solveMin()
562{
563 return solver(Minimum);
564}
565
566/*!
567 \internal
568 Maximize the original objective.
569*/
570qreal QSimplex::solveMax()
571{
572 return solver(Maximum);
573}
574
575/*!
576 \internal
577
578 Reads results from the simplified matrix and saves them in the
579 "result" member of each QSimplexVariable.
580*/
581void QSimplex::collectResults()
582{
583 // All variables are zero unless overridden below.
584
585 // ### Is this really needed? Is there any chance that an
586 // important variable remains as non-basic at the end of simplex?
587 for (int i = 0; i < variables.size(); ++i)
588 variables[i]->result = 0;
589
590 // Basic variables
591 // Update the variable indicated in the first column with the value
592 // in the last column.
593 for (int i = 1; i < rows; ++i) {
594 int index = valueAt(i, 0) - 1;
595 if (index < variables.size())
596 variables[index]->result = valueAt(i, columns - 1);
597 }
598}
599
600/*!
601 \internal
602
603 Looks for single-valued variables and remove them from the constraints list.
604*/
605bool QSimplex::simplifyConstraints(QList<QSimplexConstraint *> *constraints)
606{
607 QHash<QSimplexVariable *, qreal> results; // List of single-valued variables
608 bool modified = true; // Any chance more optimization exists?
609
610 while (modified) {
611 modified = false;
612
613 // For all constraints
614 QList<QSimplexConstraint *>::iterator iter = constraints->begin();
615 while (iter != constraints->end()) {
616 QSimplexConstraint *c = *iter;
617 if ((c->ratio == QSimplexConstraint::Equal) && (c->variables.count() == 1)) {
618 // Check whether this is a constraint of type Var == K
619 // If so, save its value to "results".
620 QSimplexVariable *variable = c->variables.constBegin().key();
621 qreal result = c->constant / c->variables.value(variable);
622
623 results.insert(variable, result);
624 variable->result = result;
625 variable->index = -1;
626 modified = true;
627
628 }
629
630 // Replace known values among their variables
631 QHash<QSimplexVariable *, qreal>::const_iterator r;
632 for (r = results.constBegin(); r != results.constEnd(); ++r) {
633 if (c->variables.contains(r.key())) {
634 c->constant -= r.value() * c->variables.take(r.key());
635 modified = true;
636 }
637 }
638
639 // Keep it normalized
640 if (c->constant < 0)
641 c->invert();
642
643 if (c->variables.isEmpty()) {
644 // If constraint became empty due to substitution, delete it.
645 if (c->isSatisfied() == false)
646 // We must ensure that the constraint soon to be deleted would not
647 // make the problem unfeasible if left behind. If that's the case,
648 // we return false so the simplex solver can properly report that.
649 return false;
650
651 delete c;
652 iter = constraints->erase(iter);
653 } else {
654 ++iter;
655 }
656 }
657 }
658
659 return true;
660}
661
662void QSimplexConstraint::invert()
663{
664 constant = -constant;
665 ratio = Ratio(2 - ratio);
666
667 QHash<QSimplexVariable *, qreal>::iterator iter;
668 for (iter = variables.begin(); iter != variables.end(); ++iter) {
669 iter.value() = -iter.value();
670 }
671}
672
673QT_END_NAMESPACE
674