1 /* 2 * Licensed to the Apache Software Foundation (ASF) under one or more 3 * contributor license agreements. See the NOTICE file distributed with 4 * this work for additional information regarding copyright ownership. 5 * The ASF licenses this file to You under the Apache License, Version 2.0 6 * (the "License"); you may not use this file except in compliance with 7 * the License. You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. 16 */ 17 18 package org.apache.commons.math3.optimization.linear; 19 20 import java.io.IOException; 21 import java.io.ObjectInputStream; 22 import java.io.ObjectOutputStream; 23 import java.io.Serializable; 24 import java.util.ArrayList; 25 import java.util.Collection; 26 import java.util.HashSet; 27 import java.util.List; 28 import java.util.Set; 29 import java.util.TreeSet; 30 31 import org.apache.commons.math3.linear.Array2DRowRealMatrix; 32 import org.apache.commons.math3.linear.MatrixUtils; 33 import org.apache.commons.math3.linear.RealMatrix; 34 import org.apache.commons.math3.linear.RealVector; 35 import org.apache.commons.math3.optimization.GoalType; 36 import org.apache.commons.math3.optimization.PointValuePair; 37 import org.apache.commons.math3.util.FastMath; 38 import org.apache.commons.math3.util.Precision; 39 40 /** 41 * A tableau for use in the Simplex method. 42 * 43 * <p> 44 * Example: 45 * <pre> 46 * W | Z | x1 | x2 | x- | s1 | s2 | a1 | RHS 47 * --------------------------------------------------- 48 * -1 0 0 0 0 0 0 1 0 <= phase 1 objective 49 * 0 1 -15 -10 0 0 0 0 0 <= phase 2 objective 50 * 0 0 1 0 0 1 0 0 2 <= constraint 1 51 * 0 0 0 1 0 0 1 0 3 <= constraint 2 52 * 0 0 1 1 0 0 0 1 4 <= constraint 3 53 * </pre> 54 * W: Phase 1 objective function</br> 55 * Z: Phase 2 objective function</br> 56 * x1 & x2: Decision variables</br> 57 * x-: Extra decision variable to allow for negative values</br> 58 * s1 & s2: Slack/Surplus variables</br> 59 * a1: Artificial variable</br> 60 * RHS: Right hand side</br> 61 * </p> 62 * @deprecated As of 3.1 (to be removed in 4.0). 63 * @since 2.0 64 */ 65 @Deprecated 66 class SimplexTableau implements Serializable { 67 68 /** Column label for negative vars. */ 69 private static final String NEGATIVE_VAR_COLUMN_LABEL = "x-"; 70 71 /** Default amount of error to accept in floating point comparisons (as ulps). */ 72 private static final int DEFAULT_ULPS = 10; 73 74 /** The cut-off threshold to zero-out entries. */ 75 private static final double CUTOFF_THRESHOLD = 1e-12; 76 77 /** Serializable version identifier. */ 78 private static final long serialVersionUID = -1369660067587938365L; 79 80 /** Linear objective function. */ 81 private final LinearObjectiveFunction f; 82 83 /** Linear constraints. */ 84 private final List<LinearConstraint> constraints; 85 86 /** Whether to restrict the variables to non-negative values. */ 87 private final boolean restrictToNonNegative; 88 89 /** The variables each column represents */ 90 private final List<String> columnLabels = new ArrayList<String>(); 91 92 /** Simple tableau. */ 93 private transient RealMatrix tableau; 94 95 /** Number of decision variables. */ 96 private final int numDecisionVariables; 97 98 /** Number of slack variables. */ 99 private final int numSlackVariables; 100 101 /** Number of artificial variables. */ 102 private int numArtificialVariables; 103 104 /** Amount of error to accept when checking for optimality. */ 105 private final double epsilon; 106 107 /** Amount of error to accept in floating point comparisons. */ 108 private final int maxUlps; 109 110 /** 111 * Build a tableau for a linear problem. 112 * @param f linear objective function 113 * @param constraints linear constraints 114 * @param goalType type of optimization goal: either {@link GoalType#MAXIMIZE} or {@link GoalType#MINIMIZE} 115 * @param restrictToNonNegative whether to restrict the variables to non-negative values 116 * @param epsilon amount of error to accept when checking for optimality 117 */ SimplexTableau(final LinearObjectiveFunction f, final Collection<LinearConstraint> constraints, final GoalType goalType, final boolean restrictToNonNegative, final double epsilon)118 SimplexTableau(final LinearObjectiveFunction f, 119 final Collection<LinearConstraint> constraints, 120 final GoalType goalType, final boolean restrictToNonNegative, 121 final double epsilon) { 122 this(f, constraints, goalType, restrictToNonNegative, epsilon, DEFAULT_ULPS); 123 } 124 125 /** 126 * Build a tableau for a linear problem. 127 * @param f linear objective function 128 * @param constraints linear constraints 129 * @param goalType type of optimization goal: either {@link GoalType#MAXIMIZE} or {@link GoalType#MINIMIZE} 130 * @param restrictToNonNegative whether to restrict the variables to non-negative values 131 * @param epsilon amount of error to accept when checking for optimality 132 * @param maxUlps amount of error to accept in floating point comparisons 133 */ SimplexTableau(final LinearObjectiveFunction f, final Collection<LinearConstraint> constraints, final GoalType goalType, final boolean restrictToNonNegative, final double epsilon, final int maxUlps)134 SimplexTableau(final LinearObjectiveFunction f, 135 final Collection<LinearConstraint> constraints, 136 final GoalType goalType, final boolean restrictToNonNegative, 137 final double epsilon, 138 final int maxUlps) { 139 this.f = f; 140 this.constraints = normalizeConstraints(constraints); 141 this.restrictToNonNegative = restrictToNonNegative; 142 this.epsilon = epsilon; 143 this.maxUlps = maxUlps; 144 this.numDecisionVariables = f.getCoefficients().getDimension() + 145 (restrictToNonNegative ? 0 : 1); 146 this.numSlackVariables = getConstraintTypeCounts(Relationship.LEQ) + 147 getConstraintTypeCounts(Relationship.GEQ); 148 this.numArtificialVariables = getConstraintTypeCounts(Relationship.EQ) + 149 getConstraintTypeCounts(Relationship.GEQ); 150 this.tableau = createTableau(goalType == GoalType.MAXIMIZE); 151 initializeColumnLabels(); 152 } 153 154 /** 155 * Initialize the labels for the columns. 156 */ initializeColumnLabels()157 protected void initializeColumnLabels() { 158 if (getNumObjectiveFunctions() == 2) { 159 columnLabels.add("W"); 160 } 161 columnLabels.add("Z"); 162 for (int i = 0; i < getOriginalNumDecisionVariables(); i++) { 163 columnLabels.add("x" + i); 164 } 165 if (!restrictToNonNegative) { 166 columnLabels.add(NEGATIVE_VAR_COLUMN_LABEL); 167 } 168 for (int i = 0; i < getNumSlackVariables(); i++) { 169 columnLabels.add("s" + i); 170 } 171 for (int i = 0; i < getNumArtificialVariables(); i++) { 172 columnLabels.add("a" + i); 173 } 174 columnLabels.add("RHS"); 175 } 176 177 /** 178 * Create the tableau by itself. 179 * @param maximize if true, goal is to maximize the objective function 180 * @return created tableau 181 */ createTableau(final boolean maximize)182 protected RealMatrix createTableau(final boolean maximize) { 183 184 // create a matrix of the correct size 185 int width = numDecisionVariables + numSlackVariables + 186 numArtificialVariables + getNumObjectiveFunctions() + 1; // + 1 is for RHS 187 int height = constraints.size() + getNumObjectiveFunctions(); 188 Array2DRowRealMatrix matrix = new Array2DRowRealMatrix(height, width); 189 190 // initialize the objective function rows 191 if (getNumObjectiveFunctions() == 2) { 192 matrix.setEntry(0, 0, -1); 193 } 194 int zIndex = (getNumObjectiveFunctions() == 1) ? 0 : 1; 195 matrix.setEntry(zIndex, zIndex, maximize ? 1 : -1); 196 RealVector objectiveCoefficients = 197 maximize ? f.getCoefficients().mapMultiply(-1) : f.getCoefficients(); 198 copyArray(objectiveCoefficients.toArray(), matrix.getDataRef()[zIndex]); 199 matrix.setEntry(zIndex, width - 1, 200 maximize ? f.getConstantTerm() : -1 * f.getConstantTerm()); 201 202 if (!restrictToNonNegative) { 203 matrix.setEntry(zIndex, getSlackVariableOffset() - 1, 204 getInvertedCoefficientSum(objectiveCoefficients)); 205 } 206 207 // initialize the constraint rows 208 int slackVar = 0; 209 int artificialVar = 0; 210 for (int i = 0; i < constraints.size(); i++) { 211 LinearConstraint constraint = constraints.get(i); 212 int row = getNumObjectiveFunctions() + i; 213 214 // decision variable coefficients 215 copyArray(constraint.getCoefficients().toArray(), matrix.getDataRef()[row]); 216 217 // x- 218 if (!restrictToNonNegative) { 219 matrix.setEntry(row, getSlackVariableOffset() - 1, 220 getInvertedCoefficientSum(constraint.getCoefficients())); 221 } 222 223 // RHS 224 matrix.setEntry(row, width - 1, constraint.getValue()); 225 226 // slack variables 227 if (constraint.getRelationship() == Relationship.LEQ) { 228 matrix.setEntry(row, getSlackVariableOffset() + slackVar++, 1); // slack 229 } else if (constraint.getRelationship() == Relationship.GEQ) { 230 matrix.setEntry(row, getSlackVariableOffset() + slackVar++, -1); // excess 231 } 232 233 // artificial variables 234 if ((constraint.getRelationship() == Relationship.EQ) || 235 (constraint.getRelationship() == Relationship.GEQ)) { 236 matrix.setEntry(0, getArtificialVariableOffset() + artificialVar, 1); 237 matrix.setEntry(row, getArtificialVariableOffset() + artificialVar++, 1); 238 matrix.setRowVector(0, matrix.getRowVector(0).subtract(matrix.getRowVector(row))); 239 } 240 } 241 242 return matrix; 243 } 244 245 /** 246 * Get new versions of the constraints which have positive right hand sides. 247 * @param originalConstraints original (not normalized) constraints 248 * @return new versions of the constraints 249 */ normalizeConstraints(Collection<LinearConstraint> originalConstraints)250 public List<LinearConstraint> normalizeConstraints(Collection<LinearConstraint> originalConstraints) { 251 List<LinearConstraint> normalized = new ArrayList<LinearConstraint>(originalConstraints.size()); 252 for (LinearConstraint constraint : originalConstraints) { 253 normalized.add(normalize(constraint)); 254 } 255 return normalized; 256 } 257 258 /** 259 * Get a new equation equivalent to this one with a positive right hand side. 260 * @param constraint reference constraint 261 * @return new equation 262 */ normalize(final LinearConstraint constraint)263 private LinearConstraint normalize(final LinearConstraint constraint) { 264 if (constraint.getValue() < 0) { 265 return new LinearConstraint(constraint.getCoefficients().mapMultiply(-1), 266 constraint.getRelationship().oppositeRelationship(), 267 -1 * constraint.getValue()); 268 } 269 return new LinearConstraint(constraint.getCoefficients(), 270 constraint.getRelationship(), constraint.getValue()); 271 } 272 273 /** 274 * Get the number of objective functions in this tableau. 275 * @return 2 for Phase 1. 1 for Phase 2. 276 */ getNumObjectiveFunctions()277 protected final int getNumObjectiveFunctions() { 278 return this.numArtificialVariables > 0 ? 2 : 1; 279 } 280 281 /** 282 * Get a count of constraints corresponding to a specified relationship. 283 * @param relationship relationship to count 284 * @return number of constraint with the specified relationship 285 */ getConstraintTypeCounts(final Relationship relationship)286 private int getConstraintTypeCounts(final Relationship relationship) { 287 int count = 0; 288 for (final LinearConstraint constraint : constraints) { 289 if (constraint.getRelationship() == relationship) { 290 ++count; 291 } 292 } 293 return count; 294 } 295 296 /** 297 * Get the -1 times the sum of all coefficients in the given array. 298 * @param coefficients coefficients to sum 299 * @return the -1 times the sum of all coefficients in the given array. 300 */ getInvertedCoefficientSum(final RealVector coefficients)301 protected static double getInvertedCoefficientSum(final RealVector coefficients) { 302 double sum = 0; 303 for (double coefficient : coefficients.toArray()) { 304 sum -= coefficient; 305 } 306 return sum; 307 } 308 309 /** 310 * Checks whether the given column is basic. 311 * @param col index of the column to check 312 * @return the row that the variable is basic in. null if the column is not basic 313 */ getBasicRow(final int col)314 protected Integer getBasicRow(final int col) { 315 Integer row = null; 316 for (int i = 0; i < getHeight(); i++) { 317 final double entry = getEntry(i, col); 318 if (Precision.equals(entry, 1d, maxUlps) && (row == null)) { 319 row = i; 320 } else if (!Precision.equals(entry, 0d, maxUlps)) { 321 return null; 322 } 323 } 324 return row; 325 } 326 327 /** 328 * Removes the phase 1 objective function, positive cost non-artificial variables, 329 * and the non-basic artificial variables from this tableau. 330 */ dropPhase1Objective()331 protected void dropPhase1Objective() { 332 if (getNumObjectiveFunctions() == 1) { 333 return; 334 } 335 336 Set<Integer> columnsToDrop = new TreeSet<Integer>(); 337 columnsToDrop.add(0); 338 339 // positive cost non-artificial variables 340 for (int i = getNumObjectiveFunctions(); i < getArtificialVariableOffset(); i++) { 341 final double entry = tableau.getEntry(0, i); 342 if (Precision.compareTo(entry, 0d, epsilon) > 0) { 343 columnsToDrop.add(i); 344 } 345 } 346 347 // non-basic artificial variables 348 for (int i = 0; i < getNumArtificialVariables(); i++) { 349 int col = i + getArtificialVariableOffset(); 350 if (getBasicRow(col) == null) { 351 columnsToDrop.add(col); 352 } 353 } 354 355 double[][] matrix = new double[getHeight() - 1][getWidth() - columnsToDrop.size()]; 356 for (int i = 1; i < getHeight(); i++) { 357 int col = 0; 358 for (int j = 0; j < getWidth(); j++) { 359 if (!columnsToDrop.contains(j)) { 360 matrix[i - 1][col++] = tableau.getEntry(i, j); 361 } 362 } 363 } 364 365 // remove the columns in reverse order so the indices are correct 366 Integer[] drop = columnsToDrop.toArray(new Integer[columnsToDrop.size()]); 367 for (int i = drop.length - 1; i >= 0; i--) { 368 columnLabels.remove((int) drop[i]); 369 } 370 371 this.tableau = new Array2DRowRealMatrix(matrix); 372 this.numArtificialVariables = 0; 373 } 374 375 /** 376 * @param src the source array 377 * @param dest the destination array 378 */ copyArray(final double[] src, final double[] dest)379 private void copyArray(final double[] src, final double[] dest) { 380 System.arraycopy(src, 0, dest, getNumObjectiveFunctions(), src.length); 381 } 382 383 /** 384 * Returns whether the problem is at an optimal state. 385 * @return whether the model has been solved 386 */ isOptimal()387 boolean isOptimal() { 388 for (int i = getNumObjectiveFunctions(); i < getWidth() - 1; i++) { 389 final double entry = tableau.getEntry(0, i); 390 if (Precision.compareTo(entry, 0d, epsilon) < 0) { 391 return false; 392 } 393 } 394 return true; 395 } 396 397 /** 398 * Get the current solution. 399 * @return current solution 400 */ getSolution()401 protected PointValuePair getSolution() { 402 int negativeVarColumn = columnLabels.indexOf(NEGATIVE_VAR_COLUMN_LABEL); 403 Integer negativeVarBasicRow = negativeVarColumn > 0 ? getBasicRow(negativeVarColumn) : null; 404 double mostNegative = negativeVarBasicRow == null ? 0 : getEntry(negativeVarBasicRow, getRhsOffset()); 405 406 Set<Integer> basicRows = new HashSet<Integer>(); 407 double[] coefficients = new double[getOriginalNumDecisionVariables()]; 408 for (int i = 0; i < coefficients.length; i++) { 409 int colIndex = columnLabels.indexOf("x" + i); 410 if (colIndex < 0) { 411 coefficients[i] = 0; 412 continue; 413 } 414 Integer basicRow = getBasicRow(colIndex); 415 if (basicRow != null && basicRow == 0) { 416 // if the basic row is found to be the objective function row 417 // set the coefficient to 0 -> this case handles unconstrained 418 // variables that are still part of the objective function 419 coefficients[i] = 0; 420 } else if (basicRows.contains(basicRow)) { 421 // if multiple variables can take a given value 422 // then we choose the first and set the rest equal to 0 423 coefficients[i] = 0 - (restrictToNonNegative ? 0 : mostNegative); 424 } else { 425 basicRows.add(basicRow); 426 coefficients[i] = 427 (basicRow == null ? 0 : getEntry(basicRow, getRhsOffset())) - 428 (restrictToNonNegative ? 0 : mostNegative); 429 } 430 } 431 return new PointValuePair(coefficients, f.getValue(coefficients)); 432 } 433 434 /** 435 * Subtracts a multiple of one row from another. 436 * <p> 437 * After application of this operation, the following will hold: 438 * <pre>minuendRow = minuendRow - multiple * subtrahendRow</pre> 439 * 440 * @param dividendRow index of the row 441 * @param divisor value of the divisor 442 */ divideRow(final int dividendRow, final double divisor)443 protected void divideRow(final int dividendRow, final double divisor) { 444 for (int j = 0; j < getWidth(); j++) { 445 tableau.setEntry(dividendRow, j, tableau.getEntry(dividendRow, j) / divisor); 446 } 447 } 448 449 /** 450 * Subtracts a multiple of one row from another. 451 * <p> 452 * After application of this operation, the following will hold: 453 * <pre>minuendRow = minuendRow - multiple * subtrahendRow</pre> 454 * 455 * @param minuendRow row index 456 * @param subtrahendRow row index 457 * @param multiple multiplication factor 458 */ subtractRow(final int minuendRow, final int subtrahendRow, final double multiple)459 protected void subtractRow(final int minuendRow, final int subtrahendRow, 460 final double multiple) { 461 for (int i = 0; i < getWidth(); i++) { 462 double result = tableau.getEntry(minuendRow, i) - tableau.getEntry(subtrahendRow, i) * multiple; 463 // cut-off values smaller than the CUTOFF_THRESHOLD, otherwise may lead to numerical instabilities 464 if (FastMath.abs(result) < CUTOFF_THRESHOLD) { 465 result = 0.0; 466 } 467 tableau.setEntry(minuendRow, i, result); 468 } 469 } 470 471 /** 472 * Get the width of the tableau. 473 * @return width of the tableau 474 */ getWidth()475 protected final int getWidth() { 476 return tableau.getColumnDimension(); 477 } 478 479 /** 480 * Get the height of the tableau. 481 * @return height of the tableau 482 */ getHeight()483 protected final int getHeight() { 484 return tableau.getRowDimension(); 485 } 486 487 /** 488 * Get an entry of the tableau. 489 * @param row row index 490 * @param column column index 491 * @return entry at (row, column) 492 */ getEntry(final int row, final int column)493 protected final double getEntry(final int row, final int column) { 494 return tableau.getEntry(row, column); 495 } 496 497 /** 498 * Set an entry of the tableau. 499 * @param row row index 500 * @param column column index 501 * @param value for the entry 502 */ setEntry(final int row, final int column, final double value)503 protected final void setEntry(final int row, final int column, 504 final double value) { 505 tableau.setEntry(row, column, value); 506 } 507 508 /** 509 * Get the offset of the first slack variable. 510 * @return offset of the first slack variable 511 */ getSlackVariableOffset()512 protected final int getSlackVariableOffset() { 513 return getNumObjectiveFunctions() + numDecisionVariables; 514 } 515 516 /** 517 * Get the offset of the first artificial variable. 518 * @return offset of the first artificial variable 519 */ getArtificialVariableOffset()520 protected final int getArtificialVariableOffset() { 521 return getNumObjectiveFunctions() + numDecisionVariables + numSlackVariables; 522 } 523 524 /** 525 * Get the offset of the right hand side. 526 * @return offset of the right hand side 527 */ getRhsOffset()528 protected final int getRhsOffset() { 529 return getWidth() - 1; 530 } 531 532 /** 533 * Get the number of decision variables. 534 * <p> 535 * If variables are not restricted to positive values, this will include 1 extra decision variable to represent 536 * the absolute value of the most negative variable. 537 * 538 * @return number of decision variables 539 * @see #getOriginalNumDecisionVariables() 540 */ getNumDecisionVariables()541 protected final int getNumDecisionVariables() { 542 return numDecisionVariables; 543 } 544 545 /** 546 * Get the original number of decision variables. 547 * @return original number of decision variables 548 * @see #getNumDecisionVariables() 549 */ getOriginalNumDecisionVariables()550 protected final int getOriginalNumDecisionVariables() { 551 return f.getCoefficients().getDimension(); 552 } 553 554 /** 555 * Get the number of slack variables. 556 * @return number of slack variables 557 */ getNumSlackVariables()558 protected final int getNumSlackVariables() { 559 return numSlackVariables; 560 } 561 562 /** 563 * Get the number of artificial variables. 564 * @return number of artificial variables 565 */ getNumArtificialVariables()566 protected final int getNumArtificialVariables() { 567 return numArtificialVariables; 568 } 569 570 /** 571 * Get the tableau data. 572 * @return tableau data 573 */ getData()574 protected final double[][] getData() { 575 return tableau.getData(); 576 } 577 578 /** {@inheritDoc} */ 579 @Override equals(Object other)580 public boolean equals(Object other) { 581 582 if (this == other) { 583 return true; 584 } 585 586 if (other instanceof SimplexTableau) { 587 SimplexTableau rhs = (SimplexTableau) other; 588 return (restrictToNonNegative == rhs.restrictToNonNegative) && 589 (numDecisionVariables == rhs.numDecisionVariables) && 590 (numSlackVariables == rhs.numSlackVariables) && 591 (numArtificialVariables == rhs.numArtificialVariables) && 592 (epsilon == rhs.epsilon) && 593 (maxUlps == rhs.maxUlps) && 594 f.equals(rhs.f) && 595 constraints.equals(rhs.constraints) && 596 tableau.equals(rhs.tableau); 597 } 598 return false; 599 } 600 601 /** {@inheritDoc} */ 602 @Override hashCode()603 public int hashCode() { 604 return Boolean.valueOf(restrictToNonNegative).hashCode() ^ 605 numDecisionVariables ^ 606 numSlackVariables ^ 607 numArtificialVariables ^ 608 Double.valueOf(epsilon).hashCode() ^ 609 maxUlps ^ 610 f.hashCode() ^ 611 constraints.hashCode() ^ 612 tableau.hashCode(); 613 } 614 615 /** 616 * Serialize the instance. 617 * @param oos stream where object should be written 618 * @throws IOException if object cannot be written to stream 619 */ writeObject(ObjectOutputStream oos)620 private void writeObject(ObjectOutputStream oos) 621 throws IOException { 622 oos.defaultWriteObject(); 623 MatrixUtils.serializeRealMatrix(tableau, oos); 624 } 625 626 /** 627 * Deserialize the instance. 628 * @param ois stream from which the object should be read 629 * @throws ClassNotFoundException if a class in the stream cannot be found 630 * @throws IOException if object cannot be read from the stream 631 */ readObject(ObjectInputStream ois)632 private void readObject(ObjectInputStream ois) 633 throws ClassNotFoundException, IOException { 634 ois.defaultReadObject(); 635 MatrixUtils.deserializeRealMatrix(this, "tableau", ois); 636 } 637 } 638