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

source code of qtbase/src/widgets/graphicsview/qsimplex_p.cpp