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 QSimplexVariable *slack;
160 QSimplexVariable *surplus;
161 QSimplexVariable *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 QSimplexVariable;
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 QSimplexVariable;
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 QSimplexVariable;
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 QHash<QSimplexVariable *, qreal>::const_iterator iter;
234 for (iter = c->variables.constBegin();
235 iter != c->variables.constEnd();
236 ++iter) {
237 setValueAt(rowIndex: i, columnIndex: iter.key()->index, value: iter.value());
238 }
239
240 setValueAt(rowIndex: i, columnIndex: columns - 1, value: c->constant);
241 }
242
243 // Set objective for the first-phase Simplex.
244 // Z = -1 * sum_of_artificial_vars
245 for (int j = firstArtificial; j < columns - 1; ++j)
246 setValueAt(rowIndex: 0, columnIndex: j, value: 1.0);
247
248 // Maximize our objective (artificial vars go to zero)
249 solveMaxHelper();
250
251 // If there is a solution where the sum of all artificial
252 // variables is zero, then all of them can be removed and yet
253 // we will have a feasible (but not optimal) solution for the
254 // original problem.
255 // Otherwise, we clean up our structures and report there is
256 // no feasible solution.
257 if ((valueAt(rowIndex: 0, columnIndex: columns - 1) != 0.0) && (qAbs(t: valueAt(rowIndex: 0, columnIndex: columns - 1)) > 0.00001)) {
258 qWarning(msg: "QSimplex: No feasible solution!");
259 clearDataStructures();
260 return false;
261 }
262
263 // Remove artificial variables. We already have a feasible
264 // solution for the first problem, thus we don't need them
265 // anymore.
266 clearColumns(first: firstArtificial, last: columns - 2);
267
268 return true;
269}
270
271/*!
272 \internal
273
274 Run simplex on the current matrix with the current objective.
275
276 This is the iterative method. The matrix lines are combined
277 as to modify the variable values towards the best solution possible.
278 The method returns when the matrix is in the optimal state.
279*/
280void QSimplex::solveMaxHelper()
281{
282 reducedRowEchelon();
283 while (iterate()) ;
284}
285
286/*!
287 \internal
288*/
289void QSimplex::setObjective(QSimplexConstraint *newObjective)
290{
291 objective = newObjective;
292}
293
294/*!
295 \internal
296*/
297void QSimplex::clearRow(int rowIndex)
298{
299 qreal *item = matrix + rowIndex * columns;
300 for (int i = 0; i < columns; ++i)
301 item[i] = 0.0;
302}
303
304/*!
305 \internal
306*/
307void QSimplex::clearColumns(int first, int last)
308{
309 for (int i = 0; i < rows; ++i) {
310 qreal *row = matrix + i * columns;
311 for (int j = first; j <= last; ++j)
312 row[j] = 0.0;
313 }
314}
315
316/*!
317 \internal
318*/
319void QSimplex::dumpMatrix()
320{
321 qDebug(msg: "---- Simplex Matrix ----\n");
322
323 QString str(" "_L1);
324 for (int j = 0; j < columns; ++j)
325 str += QString::fromLatin1(ba: " <%1 >").arg(a: j, fieldWidth: 2);
326 qDebug(msg: "%s", qPrintable(str));
327 for (int i = 0; i < rows; ++i) {
328 str = QString::fromLatin1(ba: "Row %1:").arg(a: i, fieldWidth: 2);
329
330 qreal *row = matrix + i * columns;
331 for (int j = 0; j < columns; ++j)
332 str += QString::fromLatin1(ba: "%1").arg(a: row[j], fieldWidth: 7, format: 'f', precision: 2);
333 qDebug(msg: "%s", qPrintable(str));
334 }
335 qDebug(msg: "------------------------\n");
336}
337
338/*!
339 \internal
340*/
341void QSimplex::combineRows(int toIndex, int fromIndex, qreal factor)
342{
343 if (!factor)
344 return;
345
346 qreal *from = matrix + fromIndex * columns;
347 qreal *to = matrix + toIndex * columns;
348
349 for (int j = 1; j < columns; ++j) {
350 qreal value = from[j];
351
352 // skip to[j] = to[j] + factor*0.0
353 if (value == 0.0)
354 continue;
355
356 to[j] += factor * value;
357
358 // ### Avoid Numerical errors
359 if (qAbs(t: to[j]) < 0.0000000001)
360 to[j] = 0.0;
361 }
362}
363
364/*!
365 \internal
366*/
367int QSimplex::findPivotColumn()
368{
369 qreal min = 0;
370 int minIndex = -1;
371
372 for (int j = 0; j < columns-1; ++j) {
373 if (valueAt(rowIndex: 0, columnIndex: j) < min) {
374 min = valueAt(rowIndex: 0, columnIndex: j);
375 minIndex = j;
376 }
377 }
378
379 return minIndex;
380}
381
382/*!
383 \internal
384
385 For a given pivot column, find the pivot row. That is, the row with the
386 minimum associated "quotient" where:
387
388 - quotient is the division of the value in the last column by the value
389 in the pivot column.
390 - rows with value less or equal to zero are ignored
391 - if two rows have the same quotient, lines are chosen based on the
392 highest variable index (value in the first column)
393
394 The last condition avoids a bug where artificial variables would be
395 left behind for the second-phase simplex, and with 'good'
396 constraints would be removed before it, what would lead to incorrect
397 results.
398*/
399int QSimplex::pivotRowForColumn(int column)
400{
401 qreal min = qreal(999999999999.0); // ###
402 int minIndex = -1;
403
404 for (int i = 1; i < rows; ++i) {
405 qreal divisor = valueAt(rowIndex: i, columnIndex: column);
406 if (divisor <= 0)
407 continue;
408
409 qreal quotient = valueAt(rowIndex: i, columnIndex: columns - 1) / divisor;
410 if (quotient < min) {
411 min = quotient;
412 minIndex = i;
413 } else if ((quotient == min) && (valueAt(rowIndex: i, columnIndex: 0) > valueAt(rowIndex: minIndex, columnIndex: 0))) {
414 minIndex = i;
415 }
416 }
417
418 return minIndex;
419}
420
421/*!
422 \internal
423*/
424void QSimplex::reducedRowEchelon()
425{
426 for (int i = 1; i < rows; ++i) {
427 int factorInObjectiveRow = valueAt(rowIndex: i, columnIndex: 0);
428 combineRows(toIndex: 0, fromIndex: i, factor: -1 * valueAt(rowIndex: 0, columnIndex: factorInObjectiveRow));
429 }
430}
431
432/*!
433 \internal
434
435 Does one iteration towards a better solution for the problem.
436 See 'solveMaxHelper'.
437*/
438bool QSimplex::iterate()
439{
440 // Find Pivot column
441 int pivotColumn = findPivotColumn();
442 if (pivotColumn == -1)
443 return false;
444
445 // Find Pivot row for column
446 int pivotRow = pivotRowForColumn(column: pivotColumn);
447 if (pivotRow == -1) {
448 qWarning(msg: "QSimplex: Unbounded problem!");
449 return false;
450 }
451
452 // Normalize Pivot Row
453 qreal pivot = valueAt(rowIndex: pivotRow, columnIndex: pivotColumn);
454 if (pivot != 1.0)
455 combineRows(toIndex: pivotRow, fromIndex: pivotRow, factor: (1.0 - pivot) / pivot);
456
457 // Update other rows
458 for (int row=0; row < rows; ++row) {
459 if (row == pivotRow)
460 continue;
461
462 combineRows(toIndex: row, fromIndex: pivotRow, factor: -1 * valueAt(rowIndex: row, columnIndex: pivotColumn));
463 }
464
465 // Update first column
466 setValueAt(rowIndex: pivotRow, columnIndex: 0, value: pivotColumn);
467
468 // dumpMatrix();
469 // qDebug("------------ end of iteration --------------\n");
470 return true;
471}
472
473/*!
474 \internal
475
476 Both solveMin and solveMax are interfaces to this method.
477
478 The enum SolverFactor admits 2 values: Minimum (-1) and Maximum (+1).
479
480 This method sets the original objective and runs the second phase
481 Simplex to obtain the optimal solution for the problem. As the internal
482 simplex solver is only able to _maximize_ objectives, we handle the
483 minimization case by inverting the original objective and then
484 maximizing it.
485*/
486qreal QSimplex::solver(SolverFactor factor)
487{
488 // Remove old objective
489 clearRow(rowIndex: 0);
490
491 // Set new objective in the first row of the simplex matrix
492 qreal resultOffset = 0;
493 QHash<QSimplexVariable *, qreal>::const_iterator iter;
494 for (iter = objective->variables.constBegin();
495 iter != objective->variables.constEnd();
496 ++iter) {
497
498 // Check if the variable was removed in the simplification process.
499 // If so, we save its offset to the objective function and skip adding
500 // it to the matrix.
501 if (iter.key()->index == -1) {
502 resultOffset += iter.value() * iter.key()->result;
503 continue;
504 }
505
506 setValueAt(rowIndex: 0, columnIndex: iter.key()->index, value: -1 * factor * iter.value());
507 }
508
509 solveMaxHelper();
510 collectResults();
511
512#ifdef QT_DEBUG
513 for (int i = 0; i < constraints.size(); ++i) {
514 Q_ASSERT(constraints[i]->isSatisfied());
515 }
516#endif
517
518 // Return the value calculated by the simplex plus the value of the
519 // fixed variables.
520 return (qToUnderlying(e: factor) * valueAt(rowIndex: 0, columnIndex: columns - 1)) + resultOffset;
521}
522
523/*!
524 \internal
525 Minimize the original objective.
526*/
527qreal QSimplex::solveMin()
528{
529 return solver(factor: Minimum);
530}
531
532/*!
533 \internal
534 Maximize the original objective.
535*/
536qreal QSimplex::solveMax()
537{
538 return solver(factor: Maximum);
539}
540
541/*!
542 \internal
543
544 Reads results from the simplified matrix and saves them in the
545 "result" member of each QSimplexVariable.
546*/
547void QSimplex::collectResults()
548{
549 // All variables are zero unless overridden below.
550
551 // ### Is this really needed? Is there any chance that an
552 // important variable remains as non-basic at the end of simplex?
553 for (int i = 0; i < variables.size(); ++i)
554 variables[i]->result = 0;
555
556 // Basic variables
557 // Update the variable indicated in the first column with the value
558 // in the last column.
559 for (int i = 1; i < rows; ++i) {
560 int index = valueAt(rowIndex: i, columnIndex: 0) - 1;
561 if (index < variables.size())
562 variables[index]->result = valueAt(rowIndex: i, columnIndex: columns - 1);
563 }
564}
565
566/*!
567 \internal
568
569 Looks for single-valued variables and remove them from the constraints list.
570*/
571bool QSimplex::simplifyConstraints(QList<QSimplexConstraint *> *constraints)
572{
573 QHash<QSimplexVariable *, qreal> results; // List of single-valued variables
574 bool modified = true; // Any chance more optimization exists?
575
576 while (modified) {
577 modified = false;
578
579 // For all constraints
580 QList<QSimplexConstraint *>::iterator iter = constraints->begin();
581 while (iter != constraints->end()) {
582 QSimplexConstraint *c = *iter;
583 if ((c->ratio == QSimplexConstraint::Equal) && (c->variables.size() == 1)) {
584 // Check whether this is a constraint of type Var == K
585 // If so, save its value to "results".
586 QSimplexVariable *variable = c->variables.constBegin().key();
587 qreal result = c->constant / c->variables.value(key: variable);
588
589 results.insert(key: variable, value: result);
590 variable->result = result;
591 variable->index = -1;
592 modified = true;
593
594 }
595
596 // Replace known values among their variables
597 QHash<QSimplexVariable *, qreal>::const_iterator r;
598 for (r = results.constBegin(); r != results.constEnd(); ++r) {
599 if (c->variables.contains(key: r.key())) {
600 c->constant -= r.value() * c->variables.take(key: r.key());
601 modified = true;
602 }
603 }
604
605 // Keep it normalized
606 if (c->constant < 0)
607 c->invert();
608
609 if (c->variables.isEmpty()) {
610 // If constraint became empty due to substitution, delete it.
611 if (c->isSatisfied() == false)
612 // We must ensure that the constraint soon to be deleted would not
613 // make the problem unfeasible if left behind. If that's the case,
614 // we return false so the simplex solver can properly report that.
615 return false;
616
617 delete c;
618 iter = constraints->erase(pos: iter);
619 } else {
620 ++iter;
621 }
622 }
623 }
624
625 return true;
626}
627
628void QSimplexConstraint::invert()
629{
630 constant = -constant;
631 ratio = Ratio(2 - ratio);
632
633 QHash<QSimplexVariable *, qreal>::iterator iter;
634 for (iter = variables.begin(); iter != variables.end(); ++iter) {
635 iter.value() = -iter.value();
636 }
637}
638
639QT_END_NAMESPACE
640

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