| /* |
| * Licensed to the Apache Software Foundation (ASF) under one or more |
| * contributor license agreements. See the NOTICE file distributed with |
| * this work for additional information regarding copyright ownership. |
| * The ASF licenses this file to You under the Apache License, Version 2.0 |
| * (the "License"); you may not use this file except in compliance with |
| * the License. You may obtain a copy of the License at |
| * |
| * http://www.apache.org/licenses/LICENSE-2.0 |
| * |
| * Unless required by applicable law or agreed to in writing, software |
| * distributed under the License is distributed on an "AS IS" BASIS, |
| * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| * See the License for the specific language governing permissions and |
| * limitations under the License. |
| */ |
| |
| package org.apache.commons.math.optimization.linear; |
| |
| import java.io.IOException; |
| import java.io.ObjectInputStream; |
| import java.io.ObjectOutputStream; |
| import java.io.Serializable; |
| import java.util.ArrayList; |
| import java.util.Collection; |
| import java.util.HashSet; |
| import java.util.List; |
| import java.util.Set; |
| |
| import org.apache.commons.math.linear.Array2DRowRealMatrix; |
| import org.apache.commons.math.linear.MatrixUtils; |
| import org.apache.commons.math.linear.RealMatrix; |
| import org.apache.commons.math.linear.RealVector; |
| import org.apache.commons.math.optimization.GoalType; |
| import org.apache.commons.math.optimization.RealPointValuePair; |
| import org.apache.commons.math.util.MathUtils; |
| |
| /** |
| * A tableau for use in the Simplex method. |
| * |
| * <p> |
| * Example: |
| * <pre> |
| * W | Z | x1 | x2 | x- | s1 | s2 | a1 | RHS |
| * --------------------------------------------------- |
| * -1 0 0 0 0 0 0 1 0 <= phase 1 objective |
| * 0 1 -15 -10 0 0 0 0 0 <= phase 2 objective |
| * 0 0 1 0 0 1 0 0 2 <= constraint 1 |
| * 0 0 0 1 0 0 1 0 3 <= constraint 2 |
| * 0 0 1 1 0 0 0 1 4 <= constraint 3 |
| * </pre> |
| * W: Phase 1 objective function</br> |
| * Z: Phase 2 objective function</br> |
| * x1 & x2: Decision variables</br> |
| * x-: Extra decision variable to allow for negative values</br> |
| * s1 & s2: Slack/Surplus variables</br> |
| * a1: Artificial variable</br> |
| * RHS: Right hand side</br> |
| * </p> |
| * @version $Revision: 922713 $ $Date: 2010-03-14 02:26:13 +0100 (dim. 14 mars 2010) $ |
| * @since 2.0 |
| */ |
| class SimplexTableau implements Serializable { |
| |
| /** Column label for negative vars. */ |
| private static final String NEGATIVE_VAR_COLUMN_LABEL = "x-"; |
| |
| /** Serializable version identifier. */ |
| private static final long serialVersionUID = -1369660067587938365L; |
| |
| /** Linear objective function. */ |
| private final LinearObjectiveFunction f; |
| |
| /** Linear constraints. */ |
| private final List<LinearConstraint> constraints; |
| |
| /** Whether to restrict the variables to non-negative values. */ |
| private final boolean restrictToNonNegative; |
| |
| /** The variables each column represents */ |
| private final List<String> columnLabels = new ArrayList<String>(); |
| |
| /** Simple tableau. */ |
| private transient RealMatrix tableau; |
| |
| /** Number of decision variables. */ |
| private final int numDecisionVariables; |
| |
| /** Number of slack variables. */ |
| private final int numSlackVariables; |
| |
| /** Number of artificial variables. */ |
| private int numArtificialVariables; |
| |
| /** Amount of error to accept in floating point comparisons. */ |
| private final double epsilon; |
| |
| /** |
| * Build a tableau for a linear problem. |
| * @param f linear objective function |
| * @param constraints linear constraints |
| * @param goalType type of optimization goal: either {@link GoalType#MAXIMIZE} |
| * or {@link GoalType#MINIMIZE} |
| * @param restrictToNonNegative whether to restrict the variables to non-negative values |
| * @param epsilon amount of error to accept in floating point comparisons |
| */ |
| SimplexTableau(final LinearObjectiveFunction f, |
| final Collection<LinearConstraint> constraints, |
| final GoalType goalType, final boolean restrictToNonNegative, |
| final double epsilon) { |
| this.f = f; |
| this.constraints = normalizeConstraints(constraints); |
| this.restrictToNonNegative = restrictToNonNegative; |
| this.epsilon = epsilon; |
| this.numDecisionVariables = f.getCoefficients().getDimension() + |
| (restrictToNonNegative ? 0 : 1); |
| this.numSlackVariables = getConstraintTypeCounts(Relationship.LEQ) + |
| getConstraintTypeCounts(Relationship.GEQ); |
| this.numArtificialVariables = getConstraintTypeCounts(Relationship.EQ) + |
| getConstraintTypeCounts(Relationship.GEQ); |
| this.tableau = createTableau(goalType == GoalType.MAXIMIZE); |
| initializeColumnLabels(); |
| } |
| |
| /** |
| * Initialize the labels for the columns. |
| */ |
| protected void initializeColumnLabels() { |
| if (getNumObjectiveFunctions() == 2) { |
| columnLabels.add("W"); |
| } |
| columnLabels.add("Z"); |
| for (int i = 0; i < getOriginalNumDecisionVariables(); i++) { |
| columnLabels.add("x" + i); |
| } |
| if (!restrictToNonNegative) { |
| columnLabels.add(NEGATIVE_VAR_COLUMN_LABEL); |
| } |
| for (int i = 0; i < getNumSlackVariables(); i++) { |
| columnLabels.add("s" + i); |
| } |
| for (int i = 0; i < getNumArtificialVariables(); i++) { |
| columnLabels.add("a" + i); |
| } |
| columnLabels.add("RHS"); |
| } |
| |
| /** |
| * Create the tableau by itself. |
| * @param maximize if true, goal is to maximize the objective function |
| * @return created tableau |
| */ |
| protected RealMatrix createTableau(final boolean maximize) { |
| |
| // create a matrix of the correct size |
| int width = numDecisionVariables + numSlackVariables + |
| numArtificialVariables + getNumObjectiveFunctions() + 1; // + 1 is for RHS |
| int height = constraints.size() + getNumObjectiveFunctions(); |
| Array2DRowRealMatrix matrix = new Array2DRowRealMatrix(height, width); |
| |
| // initialize the objective function rows |
| if (getNumObjectiveFunctions() == 2) { |
| matrix.setEntry(0, 0, -1); |
| } |
| int zIndex = (getNumObjectiveFunctions() == 1) ? 0 : 1; |
| matrix.setEntry(zIndex, zIndex, maximize ? 1 : -1); |
| RealVector objectiveCoefficients = |
| maximize ? f.getCoefficients().mapMultiply(-1) : f.getCoefficients(); |
| copyArray(objectiveCoefficients.getData(), matrix.getDataRef()[zIndex]); |
| matrix.setEntry(zIndex, width - 1, |
| maximize ? f.getConstantTerm() : -1 * f.getConstantTerm()); |
| |
| if (!restrictToNonNegative) { |
| matrix.setEntry(zIndex, getSlackVariableOffset() - 1, |
| getInvertedCoeffiecientSum(objectiveCoefficients)); |
| } |
| |
| // initialize the constraint rows |
| int slackVar = 0; |
| int artificialVar = 0; |
| for (int i = 0; i < constraints.size(); i++) { |
| LinearConstraint constraint = constraints.get(i); |
| int row = getNumObjectiveFunctions() + i; |
| |
| // decision variable coefficients |
| copyArray(constraint.getCoefficients().getData(), matrix.getDataRef()[row]); |
| |
| // x- |
| if (!restrictToNonNegative) { |
| matrix.setEntry(row, getSlackVariableOffset() - 1, |
| getInvertedCoeffiecientSum(constraint.getCoefficients())); |
| } |
| |
| // RHS |
| matrix.setEntry(row, width - 1, constraint.getValue()); |
| |
| // slack variables |
| if (constraint.getRelationship() == Relationship.LEQ) { |
| matrix.setEntry(row, getSlackVariableOffset() + slackVar++, 1); // slack |
| } else if (constraint.getRelationship() == Relationship.GEQ) { |
| matrix.setEntry(row, getSlackVariableOffset() + slackVar++, -1); // excess |
| } |
| |
| // artificial variables |
| if ((constraint.getRelationship() == Relationship.EQ) || |
| (constraint.getRelationship() == Relationship.GEQ)) { |
| matrix.setEntry(0, getArtificialVariableOffset() + artificialVar, 1); |
| matrix.setEntry(row, getArtificialVariableOffset() + artificialVar++, 1); |
| matrix.setRowVector(0, matrix.getRowVector(0).subtract(matrix.getRowVector(row))); |
| } |
| } |
| |
| return matrix; |
| } |
| |
| /** |
| * Get new versions of the constraints which have positive right hand sides. |
| * @param originalConstraints original (not normalized) constraints |
| * @return new versions of the constraints |
| */ |
| public List<LinearConstraint> normalizeConstraints(Collection<LinearConstraint> originalConstraints) { |
| List<LinearConstraint> normalized = new ArrayList<LinearConstraint>(); |
| for (LinearConstraint constraint : originalConstraints) { |
| normalized.add(normalize(constraint)); |
| } |
| return normalized; |
| } |
| |
| /** |
| * Get a new equation equivalent to this one with a positive right hand side. |
| * @param constraint reference constraint |
| * @return new equation |
| */ |
| private LinearConstraint normalize(final LinearConstraint constraint) { |
| if (constraint.getValue() < 0) { |
| return new LinearConstraint(constraint.getCoefficients().mapMultiply(-1), |
| constraint.getRelationship().oppositeRelationship(), |
| -1 * constraint.getValue()); |
| } |
| return new LinearConstraint(constraint.getCoefficients(), |
| constraint.getRelationship(), constraint.getValue()); |
| } |
| |
| /** |
| * Get the number of objective functions in this tableau. |
| * @return 2 for Phase 1. 1 for Phase 2. |
| */ |
| protected final int getNumObjectiveFunctions() { |
| return this.numArtificialVariables > 0 ? 2 : 1; |
| } |
| |
| /** |
| * Get a count of constraints corresponding to a specified relationship. |
| * @param relationship relationship to count |
| * @return number of constraint with the specified relationship |
| */ |
| private int getConstraintTypeCounts(final Relationship relationship) { |
| int count = 0; |
| for (final LinearConstraint constraint : constraints) { |
| if (constraint.getRelationship() == relationship) { |
| ++count; |
| } |
| } |
| return count; |
| } |
| |
| /** |
| * Get the -1 times the sum of all coefficients in the given array. |
| * @param coefficients coefficients to sum |
| * @return the -1 times the sum of all coefficients in the given array. |
| */ |
| protected static double getInvertedCoeffiecientSum(final RealVector coefficients) { |
| double sum = 0; |
| for (double coefficient : coefficients.getData()) { |
| sum -= coefficient; |
| } |
| return sum; |
| } |
| |
| /** |
| * Checks whether the given column is basic. |
| * @param col index of the column to check |
| * @return the row that the variable is basic in. null if the column is not basic |
| */ |
| protected Integer getBasicRow(final int col) { |
| Integer row = null; |
| for (int i = 0; i < getHeight(); i++) { |
| if (MathUtils.equals(getEntry(i, col), 1.0, epsilon) && (row == null)) { |
| row = i; |
| } else if (!MathUtils.equals(getEntry(i, col), 0.0, epsilon)) { |
| return null; |
| } |
| } |
| return row; |
| } |
| |
| /** |
| * Removes the phase 1 objective function, positive cost non-artificial variables, |
| * and the non-basic artificial variables from this tableau. |
| */ |
| protected void dropPhase1Objective() { |
| if (getNumObjectiveFunctions() == 1) { |
| return; |
| } |
| |
| List<Integer> columnsToDrop = new ArrayList<Integer>(); |
| columnsToDrop.add(0); |
| |
| // positive cost non-artificial variables |
| for (int i = getNumObjectiveFunctions(); i < getArtificialVariableOffset(); i++) { |
| if (MathUtils.compareTo(tableau.getEntry(0, i), 0, epsilon) > 0) { |
| columnsToDrop.add(i); |
| } |
| } |
| |
| // non-basic artificial variables |
| for (int i = 0; i < getNumArtificialVariables(); i++) { |
| int col = i + getArtificialVariableOffset(); |
| if (getBasicRow(col) == null) { |
| columnsToDrop.add(col); |
| } |
| } |
| |
| double[][] matrix = new double[getHeight() - 1][getWidth() - columnsToDrop.size()]; |
| for (int i = 1; i < getHeight(); i++) { |
| int col = 0; |
| for (int j = 0; j < getWidth(); j++) { |
| if (!columnsToDrop.contains(j)) { |
| matrix[i - 1][col++] = tableau.getEntry(i, j); |
| } |
| } |
| } |
| |
| for (int i = columnsToDrop.size() - 1; i >= 0; i--) { |
| columnLabels.remove((int) columnsToDrop.get(i)); |
| } |
| |
| this.tableau = new Array2DRowRealMatrix(matrix); |
| this.numArtificialVariables = 0; |
| } |
| |
| /** |
| * @param src the source array |
| * @param dest the destination array |
| */ |
| private void copyArray(final double[] src, final double[] dest) { |
| System.arraycopy(src, 0, dest, getNumObjectiveFunctions(), src.length); |
| } |
| |
| /** |
| * Returns whether the problem is at an optimal state. |
| * @return whether the model has been solved |
| */ |
| boolean isOptimal() { |
| for (int i = getNumObjectiveFunctions(); i < getWidth() - 1; i++) { |
| if (MathUtils.compareTo(tableau.getEntry(0, i), 0, epsilon) < 0) { |
| return false; |
| } |
| } |
| return true; |
| } |
| |
| /** |
| * Get the current solution. |
| * |
| * @return current solution |
| */ |
| protected RealPointValuePair getSolution() { |
| int negativeVarColumn = columnLabels.indexOf(NEGATIVE_VAR_COLUMN_LABEL); |
| Integer negativeVarBasicRow = negativeVarColumn > 0 ? getBasicRow(negativeVarColumn) : null; |
| double mostNegative = negativeVarBasicRow == null ? 0 : getEntry(negativeVarBasicRow, getRhsOffset()); |
| |
| Set<Integer> basicRows = new HashSet<Integer>(); |
| double[] coefficients = new double[getOriginalNumDecisionVariables()]; |
| for (int i = 0; i < coefficients.length; i++) { |
| int colIndex = columnLabels.indexOf("x" + i); |
| if (colIndex < 0) { |
| coefficients[i] = 0; |
| continue; |
| } |
| Integer basicRow = getBasicRow(colIndex); |
| if (basicRows.contains(basicRow)) { |
| // if multiple variables can take a given value |
| // then we choose the first and set the rest equal to 0 |
| coefficients[i] = 0; |
| } else { |
| basicRows.add(basicRow); |
| coefficients[i] = |
| (basicRow == null ? 0 : getEntry(basicRow, getRhsOffset())) - |
| (restrictToNonNegative ? 0 : mostNegative); |
| } |
| } |
| return new RealPointValuePair(coefficients, f.getValue(coefficients)); |
| } |
| |
| /** |
| * Subtracts a multiple of one row from another. |
| * <p> |
| * After application of this operation, the following will hold: |
| * minuendRow = minuendRow - multiple * subtrahendRow |
| * </p> |
| * @param dividendRow index of the row |
| * @param divisor value of the divisor |
| */ |
| protected void divideRow(final int dividendRow, final double divisor) { |
| for (int j = 0; j < getWidth(); j++) { |
| tableau.setEntry(dividendRow, j, tableau.getEntry(dividendRow, j) / divisor); |
| } |
| } |
| |
| /** |
| * Subtracts a multiple of one row from another. |
| * <p> |
| * After application of this operation, the following will hold: |
| * minuendRow = minuendRow - multiple * subtrahendRow |
| * </p> |
| * @param minuendRow row index |
| * @param subtrahendRow row index |
| * @param multiple multiplication factor |
| */ |
| protected void subtractRow(final int minuendRow, final int subtrahendRow, |
| final double multiple) { |
| tableau.setRowVector(minuendRow, tableau.getRowVector(minuendRow) |
| .subtract(tableau.getRowVector(subtrahendRow).mapMultiply(multiple))); |
| } |
| |
| /** |
| * Get the width of the tableau. |
| * @return width of the tableau |
| */ |
| protected final int getWidth() { |
| return tableau.getColumnDimension(); |
| } |
| |
| /** |
| * Get the height of the tableau. |
| * @return height of the tableau |
| */ |
| protected final int getHeight() { |
| return tableau.getRowDimension(); |
| } |
| |
| /** Get an entry of the tableau. |
| * @param row row index |
| * @param column column index |
| * @return entry at (row, column) |
| */ |
| protected final double getEntry(final int row, final int column) { |
| return tableau.getEntry(row, column); |
| } |
| |
| /** Set an entry of the tableau. |
| * @param row row index |
| * @param column column index |
| * @param value for the entry |
| */ |
| protected final void setEntry(final int row, final int column, |
| final double value) { |
| tableau.setEntry(row, column, value); |
| } |
| |
| /** |
| * Get the offset of the first slack variable. |
| * @return offset of the first slack variable |
| */ |
| protected final int getSlackVariableOffset() { |
| return getNumObjectiveFunctions() + numDecisionVariables; |
| } |
| |
| /** |
| * Get the offset of the first artificial variable. |
| * @return offset of the first artificial variable |
| */ |
| protected final int getArtificialVariableOffset() { |
| return getNumObjectiveFunctions() + numDecisionVariables + numSlackVariables; |
| } |
| |
| /** |
| * Get the offset of the right hand side. |
| * @return offset of the right hand side |
| */ |
| protected final int getRhsOffset() { |
| return getWidth() - 1; |
| } |
| |
| /** |
| * Get the number of decision variables. |
| * <p> |
| * If variables are not restricted to positive values, this will include 1 |
| * extra decision variable to represent the absolute value of the most |
| * negative variable. |
| * </p> |
| * @return number of decision variables |
| * @see #getOriginalNumDecisionVariables() |
| */ |
| protected final int getNumDecisionVariables() { |
| return numDecisionVariables; |
| } |
| |
| /** |
| * Get the original number of decision variables. |
| * @return original number of decision variables |
| * @see #getNumDecisionVariables() |
| */ |
| protected final int getOriginalNumDecisionVariables() { |
| return f.getCoefficients().getDimension(); |
| } |
| |
| /** |
| * Get the number of slack variables. |
| * @return number of slack variables |
| */ |
| protected final int getNumSlackVariables() { |
| return numSlackVariables; |
| } |
| |
| /** |
| * Get the number of artificial variables. |
| * @return number of artificial variables |
| */ |
| protected final int getNumArtificialVariables() { |
| return numArtificialVariables; |
| } |
| |
| /** |
| * Get the tableau data. |
| * @return tableau data |
| */ |
| protected final double[][] getData() { |
| return tableau.getData(); |
| } |
| |
| /** {@inheritDoc} */ |
| @Override |
| public boolean equals(Object other) { |
| |
| if (this == other) { |
| return true; |
| } |
| |
| if (other instanceof SimplexTableau) { |
| SimplexTableau rhs = (SimplexTableau) other; |
| return (restrictToNonNegative == rhs.restrictToNonNegative) && |
| (numDecisionVariables == rhs.numDecisionVariables) && |
| (numSlackVariables == rhs.numSlackVariables) && |
| (numArtificialVariables == rhs.numArtificialVariables) && |
| (epsilon == rhs.epsilon) && |
| f.equals(rhs.f) && |
| constraints.equals(rhs.constraints) && |
| tableau.equals(rhs.tableau); |
| } |
| return false; |
| } |
| |
| /** {@inheritDoc} */ |
| @Override |
| public int hashCode() { |
| return Boolean.valueOf(restrictToNonNegative).hashCode() ^ |
| numDecisionVariables ^ |
| numSlackVariables ^ |
| numArtificialVariables ^ |
| Double.valueOf(epsilon).hashCode() ^ |
| f.hashCode() ^ |
| constraints.hashCode() ^ |
| tableau.hashCode(); |
| } |
| |
| /** Serialize the instance. |
| * @param oos stream where object should be written |
| * @throws IOException if object cannot be written to stream |
| */ |
| private void writeObject(ObjectOutputStream oos) |
| throws IOException { |
| oos.defaultWriteObject(); |
| MatrixUtils.serializeRealMatrix(tableau, oos); |
| } |
| |
| /** Deserialize the instance. |
| * @param ois stream from which the object should be read |
| * @throws ClassNotFoundException if a class in the stream cannot be found |
| * @throws IOException if object cannot be read from the stream |
| */ |
| private void readObject(ObjectInputStream ois) |
| throws ClassNotFoundException, IOException { |
| ois.defaultReadObject(); |
| MatrixUtils.deserializeRealMatrix(this, "tableau", ois); |
| } |
| |
| } |