View Javadoc

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.math.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  
30  import org.apache.commons.math.linear.MatrixUtils;
31  import org.apache.commons.math.linear.RealMatrix;
32  import org.apache.commons.math.linear.Array2DRowRealMatrix;
33  import org.apache.commons.math.linear.RealVector;
34  import org.apache.commons.math.optimization.GoalType;
35  import org.apache.commons.math.optimization.RealPointValuePair;
36  import org.apache.commons.math.util.MathUtils;
37  
38  /**
39   * A tableau for use in the Simplex method.
40   * 
41   * <p>
42   * Example:
43   * <pre>
44   *   W |  Z |  x1 |  x2 |  x- | s1 |  s2 |  a1 |  RHS
45   * ---------------------------------------------------
46   *  -1    0    0     0     0     0     0     1     0   &lt;= phase 1 objective
47   *   0    1   -15   -10    0     0     0     0     0   &lt;= phase 2 objective
48   *   0    0    1     0     0     1     0     0     2   &lt;= constraint 1
49   *   0    0    0     1     0     0     1     0     3   &lt;= constraint 2
50   *   0    0    1     1     0     0     0     1     4   &lt;= constraint 3
51   * </pre>
52   * W: Phase 1 objective function</br>
53   * Z: Phase 2 objective function</br>
54   * x1 &amp; x2: Decision variables</br>
55   * x-: Extra decision variable to allow for negative values</br>
56   * s1 &amp; s2: Slack/Surplus variables</br>
57   * a1: Artificial variable</br>
58   * RHS: Right hand side</br>
59   * </p>
60   * @version $Revision: 797800 $ $Date: 2009-07-25 13:21:28 -0400 (Sat, 25 Jul 2009) $
61   * @since 2.0
62   */
63  class SimplexTableau implements Serializable {
64  
65      /** Serializable version identifier. */
66      private static final long serialVersionUID = -1369660067587938365L;
67  
68      /** Linear objective function. */
69      private final LinearObjectiveFunction f;
70  
71      /** Linear constraints. */
72      private final Collection<LinearConstraint> constraints;
73  
74      /** Whether to restrict the variables to non-negative values. */
75      private final boolean restrictToNonNegative;
76  
77      /** Simple tableau. */
78      protected transient RealMatrix tableau;
79  
80      /** Number of decision variables. */
81      protected final int numDecisionVariables;
82  
83      /** Number of slack variables. */
84      protected final int numSlackVariables;
85  
86      /** Number of artificial variables. */
87      protected int numArtificialVariables;
88  
89      /** Amount of error to accept in floating point comparisons. */ 
90      protected final double epsilon;
91      
92      /**
93       * Build a tableau for a linear problem.
94       * @param f linear objective function
95       * @param constraints linear constraints
96       * @param goalType type of optimization goal: either {@link GoalType#MAXIMIZE}
97       * or {@link GoalType#MINIMIZE}
98       * @param restrictToNonNegative whether to restrict the variables to non-negative values
99       * @param epsilon amount of error to accept in floating point comparisons
100      */
101     SimplexTableau(final LinearObjectiveFunction f,
102                    final Collection<LinearConstraint> constraints,
103                    final GoalType goalType, final boolean restrictToNonNegative,
104                    final double epsilon) {
105         this.f                      = f;
106         this.constraints            = constraints;
107         this.restrictToNonNegative  = restrictToNonNegative;
108         this.epsilon                = epsilon;
109         this.numDecisionVariables   = getNumVariables() + (restrictToNonNegative ? 0 : 1);
110         this.numSlackVariables      = getConstraintTypeCounts(Relationship.LEQ) +
111                                       getConstraintTypeCounts(Relationship.GEQ);
112         this.numArtificialVariables = getConstraintTypeCounts(Relationship.EQ) +
113                                       getConstraintTypeCounts(Relationship.GEQ);
114         this.tableau = new Array2DRowRealMatrix(createTableau(goalType == GoalType.MAXIMIZE));
115         initialize();
116     }
117 
118     /**
119      * Create the tableau by itself.
120      * @param maximize if true, goal is to maximize the objective function
121      * @return created tableau
122      */
123     protected double[][] createTableau(final boolean maximize) {
124 
125         // create a matrix of the correct size
126         List<LinearConstraint> constraints = getNormalizedConstraints();
127         int width = numDecisionVariables + numSlackVariables +
128         numArtificialVariables + getNumObjectiveFunctions() + 1; // + 1 is for RHS
129         int height = constraints.size() + getNumObjectiveFunctions();
130         double[][] matrix = new double[height][width];
131 
132         // initialize the objective function rows
133         if (getNumObjectiveFunctions() == 2) {
134             matrix[0][0] = -1;
135         }
136         int zIndex = (getNumObjectiveFunctions() == 1) ? 0 : 1;
137         matrix[zIndex][zIndex] = maximize ? 1 : -1;
138         RealVector objectiveCoefficients =
139             maximize ? f.getCoefficients().mapMultiply(-1) : f.getCoefficients();
140             copyArray(objectiveCoefficients.getData(), matrix[zIndex], getNumObjectiveFunctions());
141             matrix[zIndex][width - 1] =
142                 maximize ? f.getConstantTerm() : -1 * f.getConstantTerm();
143 
144                 if (!restrictToNonNegative) {
145                     matrix[zIndex][getSlackVariableOffset() - 1] =
146                         getInvertedCoeffiecientSum(objectiveCoefficients);
147                 }
148 
149                 // initialize the constraint rows
150                 int slackVar = 0;
151                 int artificialVar = 0;
152                 for (int i = 0; i < constraints.size(); i++) {
153                     LinearConstraint constraint = constraints.get(i);
154                     int row = getNumObjectiveFunctions() + i;
155 
156                     // decision variable coefficients
157                     copyArray(constraint.getCoefficients().getData(), matrix[row], 1);
158 
159                     // x-
160                     if (!restrictToNonNegative) {
161                         matrix[row][getSlackVariableOffset() - 1] =
162                             getInvertedCoeffiecientSum(constraint.getCoefficients());
163                     }
164 
165                     // RHS
166                     matrix[row][width - 1] = constraint.getValue();
167 
168                     // slack variables
169                     if (constraint.getRelationship() == Relationship.LEQ) {
170                         matrix[row][getSlackVariableOffset() + slackVar++] = 1;  // slack
171                     } else if (constraint.getRelationship() == Relationship.GEQ) {
172                         matrix[row][getSlackVariableOffset() + slackVar++] = -1; // excess
173                     }
174 
175                     // artificial variables
176                     if ((constraint.getRelationship() == Relationship.EQ) ||
177                         (constraint.getRelationship() == Relationship.GEQ)) {
178                         matrix[0][getArtificialVariableOffset() + artificialVar] = 1; 
179                         matrix[row][getArtificialVariableOffset() + artificialVar++] = 1; 
180                     }
181                 }
182 
183                 return matrix;
184     }
185 
186     /** Get the number of variables.
187      * @return number of variables
188      */
189     public int getNumVariables() {
190         return f.getCoefficients().getDimension();
191     }
192 
193     /**
194      * Get new versions of the constraints which have positive right hand sides.
195      * @return new versions of the constraints
196      */
197     public List<LinearConstraint> getNormalizedConstraints() {
198         List<LinearConstraint> normalized = new ArrayList<LinearConstraint>();
199         for (LinearConstraint constraint : constraints) {
200             normalized.add(normalize(constraint));
201         }
202         return normalized;
203     }
204 
205     /**
206      * Get a new equation equivalent to this one with a positive right hand side.
207      * @param constraint reference constraint
208      * @return new equation
209      */
210     private LinearConstraint normalize(final LinearConstraint constraint) {
211         if (constraint.getValue() < 0) {
212             return new LinearConstraint(constraint.getCoefficients().mapMultiply(-1),
213                                         constraint.getRelationship().oppositeRelationship(),
214                                         -1 * constraint.getValue());
215         }
216         return new LinearConstraint(constraint.getCoefficients(), 
217                                     constraint.getRelationship(), constraint.getValue());
218     }
219 
220     /**
221      * Get the number of objective functions in this tableau.
222      * @return 2 for Phase 1.  1 for Phase 2.
223      */
224     protected final int getNumObjectiveFunctions() {
225         return this.numArtificialVariables > 0 ? 2 : 1;
226     }
227 
228     /**
229      * Get a count of constraints corresponding to a specified relationship.
230      * @param relationship relationship to count
231      * @return number of constraint with the specified relationship
232      */
233     private int getConstraintTypeCounts(final Relationship relationship) {
234         int count = 0;
235         for (final LinearConstraint constraint : constraints) {
236             if (constraint.getRelationship() == relationship) {
237                 ++count;
238             }
239         }
240         return count;
241     }
242 
243     /**
244      * Puts the tableau in proper form by zeroing out the artificial variables
245      * in the objective function via elementary row operations.
246      */
247     private void initialize() {
248         for (int artificialVar = 0; artificialVar < numArtificialVariables; artificialVar++) {
249             int row = getBasicRow(getArtificialVariableOffset() + artificialVar);
250             subtractRow(0, row, 1.0);
251         }
252     }
253 
254     /**
255      * Get the -1 times the sum of all coefficients in the given array.
256      * @param coefficients coefficients to sum
257      * @return the -1 times the sum of all coefficients in the given array.
258      */
259     protected static double getInvertedCoeffiecientSum(final RealVector coefficients) {
260         double sum = 0;
261         for (double coefficient : coefficients.getData()) {
262             sum -= coefficient;
263         }
264         return sum;
265     }
266 
267     /**
268      * Checks whether the given column is basic.
269      * @param col index of the column to check
270      * @return the row that the variable is basic in.  null if the column is not basic
271      */
272     private Integer getBasicRow(final int col) {
273         Integer row = null;
274         for (int i = getNumObjectiveFunctions(); i < getHeight(); i++) {
275             if (MathUtils.equals(getEntry(i, col), 1.0, epsilon) && (row == null)) {
276                 row = i;
277             } else if (!MathUtils.equals(getEntry(i, col), 0.0, epsilon)) {
278                 return null;
279             }
280         }
281         return row;
282     }
283 
284     /**
285      * Removes the phase 1 objective function and artificial variables from this tableau.
286      */
287     protected void discardArtificialVariables() {
288         if (numArtificialVariables == 0) {
289             return;
290         }
291         int width = getWidth() - numArtificialVariables - 1;
292         int height = getHeight() - 1;
293         double[][] matrix = new double[height][width];
294         for (int i = 0; i < height; i++) {
295             for (int j = 0; j < width - 1; j++) {
296                 matrix[i][j] = getEntry(i + 1, j + 1);
297             }
298             matrix[i][width - 1] = getEntry(i + 1, getRhsOffset());
299         }
300         this.tableau = new Array2DRowRealMatrix(matrix);
301         this.numArtificialVariables = 0;
302     }
303 
304 
305     /**
306      * @param src the source array
307      * @param dest the destination array
308      * @param destPos the destination position
309      */
310     private void copyArray(final double[] src, final double[] dest,
311                            final int destPos) {
312         System.arraycopy(src, 0, dest, getNumObjectiveFunctions(), src.length); 
313     }
314 
315     /**
316      * Get the current solution.
317      * 
318      * @return current solution
319      */
320     protected RealPointValuePair getSolution() {
321         double[] coefficients = new double[getOriginalNumDecisionVariables()];
322         Integer basicRow =
323             getBasicRow(getNumObjectiveFunctions() + getOriginalNumDecisionVariables());
324         double mostNegative = basicRow == null ? 0 : getEntry(basicRow, getRhsOffset());
325         Set<Integer> basicRows = new HashSet<Integer>();
326         for (int i = 0; i < coefficients.length; i++) {
327             basicRow = getBasicRow(getNumObjectiveFunctions() + i);
328             if (basicRows.contains(basicRow)) {
329                 // if multiple variables can take a given value 
330                 // then we choose the first and set the rest equal to 0
331                 coefficients[i] = 0;
332             } else {
333                 basicRows.add(basicRow);
334                 coefficients[i] =
335                     (basicRow == null ? 0 : getEntry(basicRow, getRhsOffset())) -
336                     (restrictToNonNegative ? 0 : mostNegative);
337             }
338         }
339         return new RealPointValuePair(coefficients, f.getValue(coefficients));
340     }
341 
342     /**
343      * Subtracts a multiple of one row from another.
344      * <p>
345      * After application of this operation, the following will hold:
346      *   minuendRow = minuendRow - multiple * subtrahendRow
347      * </p>
348      * @param dividendRow index of the row
349      * @param divisor value of the divisor
350      */
351     protected void divideRow(final int dividendRow, final double divisor) {
352         for (int j = 0; j < getWidth(); j++) {
353             tableau.setEntry(dividendRow, j, tableau.getEntry(dividendRow, j) / divisor);
354         }
355     }
356 
357     /**
358      * Subtracts a multiple of one row from another.
359      * <p>
360      * After application of this operation, the following will hold:
361      *   minuendRow = minuendRow - multiple * subtrahendRow
362      * </p>
363      * @param minuendRow row index
364      * @param subtrahendRow row index
365      * @param multiple multiplication factor
366      */
367     protected void subtractRow(final int minuendRow, final int subtrahendRow,
368                                final double multiple) {
369         for (int j = 0; j < getWidth(); j++) {
370             tableau.setEntry(minuendRow, j, tableau.getEntry(minuendRow, j) -
371                              multiple * tableau.getEntry(subtrahendRow, j));
372         }
373     }
374 
375     /**
376      * Get the width of the tableau.
377      * @return width of the tableau
378      */
379     protected final int getWidth() {
380         return tableau.getColumnDimension();
381     }
382 
383     /**
384      * Get the height of the tableau.
385      * @return height of the tableau
386      */
387     protected final int getHeight() {
388         return tableau.getRowDimension();
389     }
390 
391     /** Get an entry of the tableau.
392      * @param row row index
393      * @param column column index
394      * @return entry at (row, column)
395      */
396     protected final double getEntry(final int row, final int column) {
397         return tableau.getEntry(row, column);
398     }
399 
400     /** Set an entry of the tableau.
401      * @param row row index
402      * @param column column index
403      * @param value for the entry
404      */
405     protected final void setEntry(final int row, final int column,
406                                   final double value) {
407         tableau.setEntry(row, column, value);
408     }
409 
410     /**
411      * Get the offset of the first slack variable.
412      * @return offset of the first slack variable
413      */
414     protected final int getSlackVariableOffset() {
415         return getNumObjectiveFunctions() + numDecisionVariables;
416     }
417 
418     /**
419      * Get the offset of the first artificial variable.
420      * @return offset of the first artificial variable
421      */
422     protected final int getArtificialVariableOffset() {
423         return getNumObjectiveFunctions() + numDecisionVariables + numSlackVariables;
424     }
425 
426     /**
427      * Get the offset of the right hand side.
428      * @return offset of the right hand side
429      */
430     protected final int getRhsOffset() {
431         return getWidth() - 1;
432     }
433 
434     /**
435      * Get the number of decision variables.
436      * <p>
437      * If variables are not restricted to positive values, this will include 1
438      * extra decision variable to represent the absolute value of the most
439      * negative variable.
440      * </p>
441      * @return number of decision variables
442      * @see #getOriginalNumDecisionVariables()
443      */
444     protected final int getNumDecisionVariables() {
445         return numDecisionVariables;
446     }
447 
448     /**
449      * Get the original number of decision variables.
450      * @return original number of decision variables
451      * @see #getNumDecisionVariables()
452      */
453     protected final int getOriginalNumDecisionVariables() {
454         return restrictToNonNegative ? numDecisionVariables : numDecisionVariables - 1;
455     }
456 
457     /**
458      * Get the number of slack variables.
459      * @return number of slack variables
460      */
461     protected final int getNumSlackVariables() {
462         return numSlackVariables;
463     }
464 
465     /**
466      * Get the number of artificial variables.
467      * @return number of artificial variables
468      */
469     protected final int getNumArtificialVariables() {
470         return numArtificialVariables;
471     }
472 
473     /**
474      * Get the tableau data.
475      * @return tableau data
476      */
477     protected final double[][] getData() {
478         return tableau.getData();
479     }
480 
481     /** {@inheritDoc} */
482     @Override
483     public boolean equals(Object other) {
484 
485       if (this == other) { 
486         return true;
487       }
488 
489       if (other == null) {
490         return false;
491       }
492 
493       try {
494 
495           SimplexTableau rhs = (SimplexTableau) other;
496           return (restrictToNonNegative  == rhs.restrictToNonNegative) &&
497                  (numDecisionVariables   == rhs.numDecisionVariables) &&
498                  (numSlackVariables      == rhs.numSlackVariables) &&
499                  (numArtificialVariables == rhs.numArtificialVariables) &&
500                  (epsilon                == rhs.epsilon) &&
501                  f.equals(rhs.f) &&
502                  constraints.equals(rhs.constraints) &&
503                  tableau.equals(rhs.tableau);
504 
505       } catch (ClassCastException ex) {
506           // ignore exception
507           return false;
508       }
509 
510     }
511     
512     /** {@inheritDoc} */
513     @Override
514     public int hashCode() {
515         return Boolean.valueOf(restrictToNonNegative).hashCode() ^
516                numDecisionVariables ^
517                numSlackVariables ^
518                numArtificialVariables ^
519                Double.valueOf(epsilon).hashCode() ^
520                f.hashCode() ^
521                constraints.hashCode() ^
522                tableau.hashCode();
523     }
524 
525     /** Serialize the instance.
526      * @param oos stream where object should be written
527      * @throws IOException if object cannot be written to stream
528      */
529     private void writeObject(ObjectOutputStream oos)
530         throws IOException {
531         oos.defaultWriteObject();
532         MatrixUtils.serializeRealMatrix(tableau, oos);
533     }
534 
535     /** Deserialize the instance.
536      * @param ois stream from which the object should be read
537      * @throws ClassNotFoundException if a class in the stream cannot be found
538      * @throws IOException if object cannot be read from the stream
539      */
540     private void readObject(ObjectInputStream ois)
541       throws ClassNotFoundException, IOException {
542         ois.defaultReadObject();
543         MatrixUtils.deserializeRealMatrix(this, "tableau", ois);
544     }
545 }