1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22 package org.hipparchus.optim.linear;
23
24 import java.io.IOException;
25 import java.io.ObjectInputStream;
26 import java.io.ObjectOutputStream;
27 import java.io.Serializable;
28 import java.util.ArrayList;
29 import java.util.Arrays;
30 import java.util.Collection;
31 import java.util.HashSet;
32 import java.util.List;
33 import java.util.Set;
34 import java.util.TreeSet;
35
36 import org.hipparchus.exception.LocalizedCoreFormats;
37 import org.hipparchus.exception.MathIllegalArgumentException;
38 import org.hipparchus.linear.Array2DRowRealMatrix;
39 import org.hipparchus.linear.RealVector;
40 import org.hipparchus.optim.PointValuePair;
41 import org.hipparchus.optim.nonlinear.scalar.GoalType;
42 import org.hipparchus.util.Precision;
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67 class SimplexTableau implements Serializable {
68
69
70 private static final String NEGATIVE_VAR_COLUMN_LABEL = "x-";
71
72
73 private static final long serialVersionUID = -1369660067587938365L;
74
75
76 private final LinearObjectiveFunction f;
77
78
79 private final List<LinearConstraint> constraints;
80
81
82 private final boolean restrictToNonNegative;
83
84
85 private final List<String> columnLabels;
86
87
88 private transient Array2DRowRealMatrix tableau;
89
90
91 private final int numDecisionVariables;
92
93
94 private final int numSlackVariables;
95
96
97 private int numArtificialVariables;
98
99
100 private final double epsilon;
101
102
103 private final int maxUlps;
104
105
106 private int[] basicVariables;
107
108
109 private int[] basicRows;
110
111
112
113
114
115
116
117
118
119
120
121
122
123 SimplexTableau(final LinearObjectiveFunction f,
124 final Collection<LinearConstraint> constraints,
125 final GoalType goalType,
126 final boolean restrictToNonNegative,
127 final double epsilon) {
128 this(f, constraints, goalType, restrictToNonNegative, epsilon, SimplexSolver.DEFAULT_ULPS);
129 }
130
131
132
133
134
135
136
137
138
139
140
141
142 SimplexTableau(final LinearObjectiveFunction f,
143 final Collection<LinearConstraint> constraints,
144 final GoalType goalType,
145 final boolean restrictToNonNegative,
146 final double epsilon,
147 final int maxUlps) throws MathIllegalArgumentException {
148 checkDimensions(f, constraints);
149 this.f = f;
150 this.constraints = normalizeConstraints(constraints);
151 this.restrictToNonNegative = restrictToNonNegative;
152 this.columnLabels = new ArrayList<>();
153 this.epsilon = epsilon;
154 this.maxUlps = maxUlps;
155 this.numDecisionVariables = f.getCoefficients().getDimension() + (restrictToNonNegative ? 0 : 1);
156 this.numSlackVariables = getConstraintTypeCounts(Relationship.LEQ) +
157 getConstraintTypeCounts(Relationship.GEQ);
158 this.numArtificialVariables = getConstraintTypeCounts(Relationship.EQ) +
159 getConstraintTypeCounts(Relationship.GEQ);
160 this.tableau = createTableau(goalType == GoalType.MAXIMIZE);
161
162
163 initializeBasicVariables(getSlackVariableOffset());
164 initializeColumnLabels();
165 }
166
167
168
169
170
171
172
173
174 private void checkDimensions(final LinearObjectiveFunction objectiveFunction,
175 final Collection<LinearConstraint> c) {
176 final int dimension = objectiveFunction.getCoefficients().getDimension();
177 for (final LinearConstraint constraint : c) {
178 final int constraintDimension = constraint.getCoefficients().getDimension();
179 if (constraintDimension != dimension) {
180 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
181 constraintDimension, dimension);
182 }
183 }
184 }
185
186
187
188 protected void initializeColumnLabels() {
189 if (getNumObjectiveFunctions() == 2) {
190 columnLabels.add("W");
191 }
192 columnLabels.add("Z");
193 for (int i = 0; i < getOriginalNumDecisionVariables(); i++) {
194 columnLabels.add("x" + i);
195 }
196 if (!restrictToNonNegative) {
197 columnLabels.add(NEGATIVE_VAR_COLUMN_LABEL);
198 }
199 for (int i = 0; i < getNumSlackVariables(); i++) {
200 columnLabels.add("s" + i);
201 }
202 for (int i = 0; i < getNumArtificialVariables(); i++) {
203 columnLabels.add("a" + i);
204 }
205 columnLabels.add("RHS");
206 }
207
208
209
210
211
212
213 protected Array2DRowRealMatrix createTableau(final boolean maximize) {
214
215
216 int width = numDecisionVariables + numSlackVariables +
217 numArtificialVariables + getNumObjectiveFunctions() + 1;
218 int height = constraints.size() + getNumObjectiveFunctions();
219 Array2DRowRealMatrix matrix = new Array2DRowRealMatrix(height, width);
220
221
222 if (getNumObjectiveFunctions() == 2) {
223 matrix.setEntry(0, 0, -1);
224 }
225
226 int zIndex = (getNumObjectiveFunctions() == 1) ? 0 : 1;
227 matrix.setEntry(zIndex, zIndex, maximize ? 1 : -1);
228 RealVector objectiveCoefficients = maximize ? f.getCoefficients().mapMultiply(-1) : f.getCoefficients();
229 copyArray(objectiveCoefficients.toArray(), matrix.getDataRef()[zIndex]);
230 matrix.setEntry(zIndex, width - 1, maximize ? f.getConstantTerm() : -1 * f.getConstantTerm());
231
232 if (!restrictToNonNegative) {
233 matrix.setEntry(zIndex, getSlackVariableOffset() - 1,
234 getInvertedCoefficientSum(objectiveCoefficients));
235 }
236
237
238 int slackVar = 0;
239 int artificialVar = 0;
240 for (int i = 0; i < constraints.size(); i++) {
241 LinearConstraint constraint = constraints.get(i);
242 int row = getNumObjectiveFunctions() + i;
243
244
245 copyArray(constraint.getCoefficients().toArray(), matrix.getDataRef()[row]);
246
247
248 if (!restrictToNonNegative) {
249 matrix.setEntry(row, getSlackVariableOffset() - 1,
250 getInvertedCoefficientSum(constraint.getCoefficients()));
251 }
252
253
254 matrix.setEntry(row, width - 1, constraint.getValue());
255
256
257 if (constraint.getRelationship() == Relationship.LEQ) {
258 matrix.setEntry(row, getSlackVariableOffset() + slackVar++, 1);
259 } else if (constraint.getRelationship() == Relationship.GEQ) {
260 matrix.setEntry(row, getSlackVariableOffset() + slackVar++, -1);
261 }
262
263
264 if ((constraint.getRelationship() == Relationship.EQ) ||
265 (constraint.getRelationship() == Relationship.GEQ)) {
266 matrix.setEntry(0, getArtificialVariableOffset() + artificialVar, 1);
267 matrix.setEntry(row, getArtificialVariableOffset() + artificialVar++, 1);
268 matrix.setRowVector(0, matrix.getRowVector(0).subtract(matrix.getRowVector(row)));
269 }
270 }
271
272 return matrix;
273 }
274
275
276
277
278
279
280 public List<LinearConstraint> normalizeConstraints(Collection<LinearConstraint> originalConstraints) {
281 List<LinearConstraint> normalized = new ArrayList<>(originalConstraints.size());
282 for (LinearConstraint constraint : originalConstraints) {
283 normalized.add(normalize(constraint));
284 }
285 return normalized;
286 }
287
288
289
290
291
292
293 private LinearConstraint normalize(final LinearConstraint constraint) {
294 if (constraint.getValue() < 0) {
295 return new LinearConstraint(constraint.getCoefficients().mapMultiply(-1),
296 constraint.getRelationship().oppositeRelationship(),
297 -1 * constraint.getValue());
298 }
299 return new LinearConstraint(constraint.getCoefficients(),
300 constraint.getRelationship(), constraint.getValue());
301 }
302
303
304
305
306
307 protected final int getNumObjectiveFunctions() {
308 return this.numArtificialVariables > 0 ? 2 : 1;
309 }
310
311
312
313
314
315
316 private int getConstraintTypeCounts(final Relationship relationship) {
317 int count = 0;
318 for (final LinearConstraint constraint : constraints) {
319 if (constraint.getRelationship() == relationship) {
320 ++count;
321 }
322 }
323 return count;
324 }
325
326
327
328
329
330
331 protected static double getInvertedCoefficientSum(final RealVector coefficients) {
332 double sum = 0;
333 for (double coefficient : coefficients.toArray()) {
334 sum -= coefficient;
335 }
336 return sum;
337 }
338
339
340
341
342
343
344 protected Integer getBasicRow(final int col) {
345 final int row = basicVariables[col];
346 return row == -1 ? null : row;
347 }
348
349
350
351
352
353
354 protected int getBasicVariable(final int row) {
355 return basicRows[row];
356 }
357
358
359
360
361
362 private void initializeBasicVariables(final int startColumn) {
363 basicVariables = new int[getWidth() - 1];
364 basicRows = new int[getHeight()];
365
366 Arrays.fill(basicVariables, -1);
367
368 for (int i = startColumn; i < getWidth() - 1; i++) {
369 Integer row = findBasicRow(i);
370 if (row != null) {
371 basicVariables[i] = row;
372 basicRows[row] = i;
373 }
374 }
375 }
376
377
378
379
380
381
382 private Integer findBasicRow(final int col) {
383 Integer row = null;
384 for (int i = 0; i < getHeight(); i++) {
385 final double entry = getEntry(i, col);
386 if (Precision.equals(entry, 1d, maxUlps) && (row == null)) {
387 row = i;
388 } else if (!Precision.equals(entry, 0d, maxUlps)) {
389 return null;
390 }
391 }
392 return row;
393 }
394
395
396
397
398
399 protected void dropPhase1Objective() {
400 if (getNumObjectiveFunctions() == 1) {
401 return;
402 }
403
404 final Set<Integer> columnsToDrop = new TreeSet<>();
405 columnsToDrop.add(0);
406
407
408 for (int i = getNumObjectiveFunctions(); i < getArtificialVariableOffset(); i++) {
409 final double entry = getEntry(0, i);
410 if (Precision.compareTo(entry, 0d, epsilon) > 0) {
411 columnsToDrop.add(i);
412 }
413 }
414
415
416 for (int i = 0; i < getNumArtificialVariables(); i++) {
417 int col = i + getArtificialVariableOffset();
418 if (getBasicRow(col) == null) {
419 columnsToDrop.add(col);
420 }
421 }
422
423 final double[][] matrix = new double[getHeight() - 1][getWidth() - columnsToDrop.size()];
424 for (int i = 1; i < getHeight(); i++) {
425 int col = 0;
426 for (int j = 0; j < getWidth(); j++) {
427 if (!columnsToDrop.contains(j)) {
428 matrix[i - 1][col++] = getEntry(i, j);
429 }
430 }
431 }
432
433
434 Integer[] drop = columnsToDrop.toArray(new Integer[0]);
435 for (int i = drop.length - 1; i >= 0; i--) {
436 columnLabels.remove((int) drop[i]);
437 }
438
439 this.tableau = new Array2DRowRealMatrix(matrix);
440 this.numArtificialVariables = 0;
441
442 initializeBasicVariables(getNumObjectiveFunctions());
443 }
444
445
446
447
448
449 private void copyArray(final double[] src, final double[] dest) {
450 System.arraycopy(src, 0, dest, getNumObjectiveFunctions(), src.length);
451 }
452
453
454
455
456
457 boolean isOptimal() {
458 final double[] objectiveFunctionRow = getRow(0);
459 final int end = getRhsOffset();
460 for (int i = getNumObjectiveFunctions(); i < end; i++) {
461 final double entry = objectiveFunctionRow[i];
462 if (Precision.compareTo(entry, 0d, epsilon) < 0) {
463 return false;
464 }
465 }
466 return true;
467 }
468
469
470
471
472
473 protected PointValuePair getSolution() {
474 int negativeVarColumn = columnLabels.indexOf(NEGATIVE_VAR_COLUMN_LABEL);
475 Integer negativeVarBasicRow = negativeVarColumn > 0 ? getBasicRow(negativeVarColumn) : null;
476 double mostNegative = negativeVarBasicRow == null ? 0 : getEntry(negativeVarBasicRow, getRhsOffset());
477
478 final Set<Integer> usedBasicRows = new HashSet<>();
479 final double[] coefficients = new double[getOriginalNumDecisionVariables()];
480 for (int i = 0; i < coefficients.length; i++) {
481 int colIndex = columnLabels.indexOf("x" + i);
482 if (colIndex < 0) {
483 coefficients[i] = 0;
484 continue;
485 }
486 Integer basicRow = getBasicRow(colIndex);
487 if (basicRow != null && basicRow == 0) {
488
489
490
491 coefficients[i] = 0;
492 } else if (usedBasicRows.contains(basicRow)) {
493
494
495 coefficients[i] = 0 - (restrictToNonNegative ? 0 : mostNegative);
496 } else {
497 usedBasicRows.add(basicRow);
498 coefficients[i] =
499 (basicRow == null ? 0 : getEntry(basicRow, getRhsOffset())) -
500 (restrictToNonNegative ? 0 : mostNegative);
501 }
502 }
503 return new PointValuePair(coefficients, f.value(coefficients));
504 }
505
506
507
508
509
510
511
512 protected void performRowOperations(int pivotCol, int pivotRow) {
513
514 final double pivotVal = getEntry(pivotRow, pivotCol);
515 divideRow(pivotRow, pivotVal);
516
517
518 for (int i = 0; i < getHeight(); i++) {
519 if (i != pivotRow) {
520 final double multiplier = getEntry(i, pivotCol);
521 if (multiplier != 0.0) {
522 subtractRow(i, pivotRow, multiplier);
523 }
524 }
525 }
526
527
528 final int previousBasicVariable = getBasicVariable(pivotRow);
529 basicVariables[previousBasicVariable] = -1;
530 basicVariables[pivotCol] = pivotRow;
531 basicRows[pivotRow] = pivotCol;
532 }
533
534
535
536
537
538
539
540
541
542
543 protected void divideRow(final int dividendRowIndex, final double divisor) {
544 final double[] dividendRow = getRow(dividendRowIndex);
545 for (int j = 0; j < getWidth(); j++) {
546 dividendRow[j] /= divisor;
547 }
548 }
549
550
551
552
553
554
555
556
557
558
559
560 protected void subtractRow(final int minuendRowIndex, final int subtrahendRowIndex, final double multiplier) {
561 final double[] minuendRow = getRow(minuendRowIndex);
562 final double[] subtrahendRow = getRow(subtrahendRowIndex);
563 for (int i = 0; i < getWidth(); i++) {
564 minuendRow[i] -= subtrahendRow[i] * multiplier;
565 }
566 }
567
568
569
570
571
572 protected final int getWidth() {
573 return tableau.getColumnDimension();
574 }
575
576
577
578
579
580 protected final int getHeight() {
581 return tableau.getRowDimension();
582 }
583
584
585
586
587
588
589
590 protected final double getEntry(final int row, final int column) {
591 return tableau.getEntry(row, column);
592 }
593
594
595
596
597
598
599
600 protected final void setEntry(final int row, final int column, final double value) {
601 tableau.setEntry(row, column, value);
602 }
603
604
605
606
607
608 protected final int getSlackVariableOffset() {
609 return getNumObjectiveFunctions() + numDecisionVariables;
610 }
611
612
613
614
615
616 protected final int getArtificialVariableOffset() {
617 return getNumObjectiveFunctions() + numDecisionVariables + numSlackVariables;
618 }
619
620
621
622
623
624 protected final int getRhsOffset() {
625 return getWidth() - 1;
626 }
627
628
629
630
631
632
633
634
635
636
637 protected final int getNumDecisionVariables() {
638 return numDecisionVariables;
639 }
640
641
642
643
644
645
646 protected final int getOriginalNumDecisionVariables() {
647 return f.getCoefficients().getDimension();
648 }
649
650
651
652
653
654 protected final int getNumSlackVariables() {
655 return numSlackVariables;
656 }
657
658
659
660
661
662 protected final int getNumArtificialVariables() {
663 return numArtificialVariables;
664 }
665
666
667
668
669
670
671 protected final double[] getRow(int row) {
672 return tableau.getDataRef()[row];
673 }
674
675
676
677
678
679 protected final double[][] getData() {
680 return tableau.getData();
681 }
682
683
684 @Override
685 public boolean equals(Object other) {
686
687 if (this == other) {
688 return true;
689 }
690
691 if (other instanceof SimplexTableau) {
692 SimplexTableau rhs = (SimplexTableau) other;
693 return (restrictToNonNegative == rhs.restrictToNonNegative) &&
694 (numDecisionVariables == rhs.numDecisionVariables) &&
695 (numSlackVariables == rhs.numSlackVariables) &&
696 (numArtificialVariables == rhs.numArtificialVariables) &&
697 (epsilon == rhs.epsilon) &&
698 (maxUlps == rhs.maxUlps) &&
699 f.equals(rhs.f) &&
700 constraints.equals(rhs.constraints) &&
701 tableau.equals(rhs.tableau);
702 }
703 return false;
704 }
705
706
707 @Override
708 public int hashCode() {
709 return Boolean.valueOf(restrictToNonNegative).hashCode() ^
710 numDecisionVariables ^
711 numSlackVariables ^
712 numArtificialVariables ^
713 Double.valueOf(epsilon).hashCode() ^
714 maxUlps ^
715 f.hashCode() ^
716 constraints.hashCode() ^
717 tableau.hashCode();
718 }
719
720
721
722
723
724
725 private void writeObject(ObjectOutputStream oos)
726 throws IOException {
727 oos.defaultWriteObject();
728 final int n = tableau.getRowDimension();
729 final int m = tableau.getColumnDimension();
730 oos.writeInt(n);
731 oos.writeInt(m);
732 for (int i = 0; i < n; ++i) {
733 for (int j = 0; j < m; ++j) {
734 oos.writeDouble(tableau.getEntry(i, j));
735 }
736 }
737 }
738
739
740
741
742
743
744
745 private void readObject(ObjectInputStream ois)
746 throws ClassNotFoundException, IOException {
747 ois.defaultReadObject();
748
749
750 final int n = ois.readInt();
751 final int m = ois.readInt();
752 final double[][] data = new double[n][m];
753 for (int i = 0; i < n; ++i) {
754 final double[] dataI = data[i];
755 for (int j = 0; j < m; ++j) {
756 dataI[j] = ois.readDouble();
757 }
758 }
759
760
761 tableau = new Array2DRowRealMatrix(data, false);
762
763 }
764 }