added Benjamin's patch from 2009-09-07

JIRA: MATH-286

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@812390 13f79535-47bb-0310-9956-ffa450edef68

src/main/java/org/apache/commons/math/optimization/linear/SimplexSolver.java

@@ -17,6 +17,9 @@
 
 package org.apache.commons.math.optimization.linear;
 
+import java.util.ArrayList;
+import java.util.List;
+
 import org.apache.commons.math.optimization.OptimizationException;
 import org.apache.commons.math.optimization.RealPointValuePair;
 import org.apache.commons.math.util.MathUtils;
@@ -73,23 +76,42 @@ public class SimplexSolver extends AbstractLinearOptimizer {
      * @param col the column to test the ratio of.  See {@link #getPivotColumn(SimplexTableau)}
      * @return row with the minimum ratio
      */
-    private Integer getPivotRow(final int col, final SimplexTableau tableau) {
+    private Integer getPivotRow(SimplexTableau tableau, final int col) {
+        // create a list of all the rows that tie for the lowest score in the minimum ratio test
+        List<Integer> minRatioPositions = new ArrayList<Integer>();
         double minRatio = Double.MAX_VALUE;
-        Integer minRatioPos = null;
         for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getHeight(); i++) {
             final double rhs = tableau.getEntry(i, tableau.getWidth() - 1);
             final double entry = tableau.getEntry(i, col);
             if (MathUtils.compareTo(entry, 0, epsilon) > 0) {
                 final double ratio = rhs / entry;
-                if (ratio < minRatio) {
+                if (MathUtils.equals(ratio, minRatio, epsilon)) {
+                    minRatioPositions.add(i);
+                } else if (ratio < minRatio) {
                     minRatio = ratio;
-                    minRatioPos = i;
+                    minRatioPositions = new ArrayList<Integer>();
+                    minRatioPositions.add(i);
                 }
             }
         }
-        return minRatioPos;
-    }
 
+        if (minRatioPositions.size() == 0) {
+          return null;
+        } else if (minRatioPositions.size() > 1) {
+          // there's a degeneracy as indicated by a tie in the minimum ratio test
+          // check if there's an artificial variable that can be forced out of the basis
+          for (Integer row : minRatioPositions) {
+            for (int i = 0; i < tableau.getNumArtificialVariables(); i++) {
+              int column = i + tableau.getArtificialVariableOffset();
+              if (MathUtils.equals(tableau.getEntry(row, column), 1, epsilon) &&
+                  row.equals(tableau.getBasicRow(column))) {
+                return row;
+              }
+            }
+          }
+        }
+        return minRatioPositions.get(0);
+    }
 
     /**
      * Runs one iteration of the Simplex method on the given model.
@@ -103,7 +125,7 @@ public class SimplexSolver extends AbstractLinearOptimizer {
         incrementIterationsCounter();
 
         Integer pivotCol = getPivotColumn(tableau);
-        Integer pivotRow = getPivotRow(pivotCol, tableau);
+        Integer pivotRow = getPivotRow(tableau, pivotCol);
         if (pivotRow == null) {
             throw new UnboundedSolutionException();
         }
@@ -121,40 +143,6 @@ public class SimplexSolver extends AbstractLinearOptimizer {
         }
     }
 
-    /**
-     * Checks whether Phase 1 is solved.
-     * @param tableau simple tableau for the problem
-     * @return whether Phase 1 is solved
-     */
-    private boolean isPhase1Solved(final SimplexTableau tableau) {
-        if (tableau.getNumArtificialVariables() == 0) {
-            return true;
-        }
-        for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getWidth() - 1; i++) {
-            if (MathUtils.compareTo(tableau.getEntry(0, i), 0, epsilon) < 0) {
-                return false;
-            }
-        }
-        return true;
-    }
-
-    /**
-     * Returns whether the problem is at an optimal state.
-     * @param tableau simple tableau for the problem
-     * @return whether the model has been solved
-     */
-    public boolean isOptimal(final SimplexTableau tableau) {
-        if (tableau.getNumArtificialVariables() > 0) {
-            return false;
-        }
-        for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getWidth() - 1; i++) {
-            if (MathUtils.compareTo(tableau.getEntry(0, i), 0, epsilon) < 0) {
-                return false;
-            }
-        }
-        return true;
-    }
-
     /**
      * Solves Phase 1 of the Simplex method.
      * @param tableau simple tableau for the problem
@@ -162,14 +150,14 @@ public class SimplexSolver extends AbstractLinearOptimizer {
      * exceeded, or if the problem is found not to have a bounded solution, or
      * if there is no feasible solution
      */
-    protected void solvePhase1(final SimplexTableau tableau)
+    protected void solvePhase1(final SimplexTableau tableau) throws OptimizationException {
-        throws OptimizationException {
+
         // make sure we're in Phase 1
         if (tableau.getNumArtificialVariables() == 0) {
             return;
         }
 
-        while (!isPhase1Solved(tableau)) {
+        while (!tableau.isOptimal()) {
             doIteration(tableau);
         }
 
@@ -181,13 +169,14 @@ public class SimplexSolver extends AbstractLinearOptimizer {
 
     /** {@inheritDoc} */
     @Override
-    public RealPointValuePair doOptimize()
+    public RealPointValuePair doOptimize() throws OptimizationException {
-        throws OptimizationException {
         final SimplexTableau tableau =
             new SimplexTableau(function, linearConstraints, goal, nonNegative, epsilon);
+
         solvePhase1(tableau);
         tableau.discardArtificialVariables();
-        while (!isOptimal(tableau)) {
+
+        while (!tableau.isOptimal()) {
             doIteration(tableau);
         }
         return tableau.getSolution();

src/main/java/org/apache/commons/math/optimization/linear/SimplexTableau.java

@@ -27,9 +27,9 @@ 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.Array2DRowRealMatrix;
 import org.apache.commons.math.linear.RealVector;
 import org.apache.commons.math.optimization.GoalType;
 import org.apache.commons.math.optimization.RealPointValuePair;
@@ -106,7 +106,8 @@ class SimplexTableau implements Serializable {
         this.constraints            = normalizeConstraints(constraints);
         this.restrictToNonNegative  = restrictToNonNegative;
         this.epsilon                = epsilon;
-        this.numDecisionVariables   = getNumVariables() + (restrictToNonNegative ? 0 : 1);
+        this.numDecisionVariables   = f.getCoefficients().getDimension() +
+                                      (restrictToNonNegative ? 0 : 1);
         this.numSlackVariables      = getConstraintTypeCounts(Relationship.LEQ) +
                                       getConstraintTypeCounts(Relationship.GEQ);
         this.numArtificialVariables = getConstraintTypeCounts(Relationship.EQ) +
@@ -183,13 +184,6 @@ class SimplexTableau implements Serializable {
 
     }
 
-    /** Get the number of variables.
-     * @return number of variables
-     */
-    public int getNumVariables() {
-        return f.getCoefficients().getDimension();
-    }
-
     /**
      * Get new versions of the constraints which have positive right hand sides.
      * @param originalConstraints original (not normalized) constraints
@@ -270,7 +264,7 @@ class SimplexTableau implements Serializable {
      * @param col index of the column to check
      * @return the row that the variable is basic in.  null if the column is not basic
      */
-    private Integer getBasicRow(final int col) {
+    Integer getBasicRow(final int col) {
         return getBasicRow(col, true);
     }
 
@@ -323,7 +317,6 @@ class SimplexTableau implements Serializable {
         this.numArtificialVariables = 0;
     }
 
-
     /**
      * @param src the source array
      * @param dest the destination array
@@ -332,6 +325,19 @@ class SimplexTableau implements Serializable {
         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.
      *

src/test/java/org/apache/commons/math/optimization/linear/SimplexSolverTest.java

@@ -17,13 +17,11 @@
 
 package org.apache.commons.math.optimization.linear;
 
-import static org.junit.Assert.assertEquals;
+import org.junit.Assert;
 
 import java.util.ArrayList;
 import java.util.Collection;
 
-import org.apache.commons.math.linear.RealVector;
-import org.apache.commons.math.linear.ArrayRealVector;
 import org.apache.commons.math.optimization.GoalType;
 import org.apache.commons.math.optimization.OptimizationException;
 import org.apache.commons.math.optimization.RealPointValuePair;
@@ -42,20 +40,34 @@ public class SimplexSolverTest {
         SimplexSolver solver = new SimplexSolver();
         RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MINIMIZE, true);
 
-        assertEquals(0.0, solution.getPoint()[0], .0000001);
+        Assert.assertEquals(0.0, solution.getPoint()[0], .0000001);
-        assertEquals(1.0, solution.getPoint()[1], .0000001);
+        Assert.assertEquals(1.0, solution.getPoint()[1], .0000001);
-        assertEquals(1.0, solution.getPoint()[2], .0000001);
+        Assert.assertEquals(1.0, solution.getPoint()[2], .0000001);
-        assertEquals(3.0, solution.getValue(), .0000001);
+        Assert.assertEquals(3.0, solution.getValue(), .0000001);
     }
 
     @Test
     public void testMath286() throws OptimizationException {
-      LinearObjectiveFunction f = new LinearObjectiveFunction(new double[] { 0.2, 0.3 }, 0 );
+        LinearObjectiveFunction f = new LinearObjectiveFunction(new double[] { 0.2, 0.3 }, 0 );
-      Collection<LinearConstraint> constraints = new ArrayList<LinearConstraint>();
+        Collection<LinearConstraint> constraints = new ArrayList<LinearConstraint>();
-      constraints.add(new LinearConstraint(new double[] { 1, 1 }, Relationship.EQ, 23.0));
+        constraints.add(new LinearConstraint(new double[] { 1, 1 }, Relationship.EQ, 23.0));
 
-      RealPointValuePair solution = new SimplexSolver().optimize(f, constraints, GoalType.MAXIMIZE, true);
+        SimplexSolver solver = new SimplexSolver();
-      assertEquals(6.9, solution.getValue(), .0000001);
+        RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MAXIMIZE, true);
+        Assert.assertEquals(6.9, solution.getValue(), .0000001);
+    }
+
+    @Test
+    public void testDegeneracy() throws OptimizationException {
+        LinearObjectiveFunction f = new LinearObjectiveFunction(new double[] { 0.8, 0.7 }, 0 );
+        Collection<LinearConstraint> constraints = new ArrayList<LinearConstraint>();
+        constraints.add(new LinearConstraint(new double[] { 1, 1 }, Relationship.LEQ, 18.0));
+        constraints.add(new LinearConstraint(new double[] { 1, 0 }, Relationship.GEQ, 10.0));
+        constraints.add(new LinearConstraint(new double[] { 0, 1 }, Relationship.GEQ, 8.0));
+
+        SimplexSolver solver = new SimplexSolver();
+        RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MAXIMIZE, true);
+        Assert.assertEquals(13.6, solution.getValue(), .0000001);
     }
 
     @Test
@@ -70,7 +82,7 @@ public class SimplexSolverTest {
 
         SimplexSolver solver = new SimplexSolver();
         RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MAXIMIZE, true);
-        assertEquals(10.0, solution.getValue(), .0000001);
+        Assert.assertEquals(10.0, solution.getValue(), .0000001);
     }
 
     @Test
@@ -80,9 +92,9 @@ public class SimplexSolverTest {
         constraints.add(new LinearConstraint(new double[] { 2, 0 }, Relationship.GEQ, -1.0));
         SimplexSolver solver = new SimplexSolver();
         RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MINIMIZE, true);
-        assertEquals(0, solution.getValue(), .0000001);
+        Assert.assertEquals(0, solution.getValue(), .0000001);
-        assertEquals(0, solution.getPoint()[0], .0000001);
+        Assert.assertEquals(0, solution.getPoint()[0], .0000001);
-        assertEquals(0, solution.getPoint()[1], .0000001);
+        Assert.assertEquals(0, solution.getPoint()[1], .0000001);
     }
 
     @Test(expected=NoFeasibleSolutionException.class)
@@ -105,9 +117,9 @@ public class SimplexSolverTest {
 
         SimplexSolver solver = new SimplexSolver();
         RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MAXIMIZE, false);
-        assertEquals(2.0, solution.getPoint()[0], 0.0);
+        Assert.assertEquals(2.0, solution.getPoint()[0], 0.0);
-        assertEquals(2.0, solution.getPoint()[1], 0.0);
+        Assert.assertEquals(2.0, solution.getPoint()[1], 0.0);
-        assertEquals(57.0, solution.getValue(), 0.0);
+        Assert.assertEquals(57.0, solution.getValue(), 0.0);
     }
 
     @Test
@@ -118,8 +130,8 @@ public class SimplexSolverTest {
 
         SimplexSolver solver = new SimplexSolver();
         RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MAXIMIZE, false);
-        assertEquals(10.0, solution.getPoint()[0], 0.0);
+        Assert.assertEquals(10.0, solution.getPoint()[0], 0.0);
-        assertEquals(30.0, solution.getValue(), 0.0);
+        Assert.assertEquals(30.0, solution.getValue(), 0.0);
     }
 
     /**
@@ -136,9 +148,9 @@ public class SimplexSolverTest {
 
         SimplexSolver solver = new SimplexSolver();
         RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MAXIMIZE, false);
-        assertEquals(2.0, solution.getPoint()[0], 0.0);
+        Assert.assertEquals(2.0, solution.getPoint()[0], 0.0);
-        assertEquals(2.0, solution.getPoint()[1], 0.0);
+        Assert.assertEquals(2.0, solution.getPoint()[1], 0.0);
-        assertEquals(50.0, solution.getValue(), 0.0);
+        Assert.assertEquals(50.0, solution.getValue(), 0.0);
     }
 
     @Test
@@ -151,9 +163,9 @@ public class SimplexSolverTest {
 
         SimplexSolver solver = new SimplexSolver();
         RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MINIMIZE, false);
-        assertEquals(4.0, solution.getPoint()[0], 0.0);
+        Assert.assertEquals(4.0, solution.getPoint()[0], 0.0);
-        assertEquals(0.0, solution.getPoint()[1], 0.0);
+        Assert.assertEquals(0.0, solution.getPoint()[1], 0.0);
-        assertEquals(-13.0, solution.getValue(), 0.0);
+        Assert.assertEquals(-13.0, solution.getValue(), 0.0);
     }
 
     @Test
@@ -165,9 +177,9 @@ public class SimplexSolverTest {
 
         SimplexSolver solver = new SimplexSolver();
         RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MAXIMIZE, false);
-        assertEquals(-2.0, solution.getPoint()[0], 0.0);
+        Assert.assertEquals(-2.0, solution.getPoint()[0], 0.0);
-        assertEquals(8.0, solution.getPoint()[1], 0.0);
+        Assert.assertEquals(8.0, solution.getPoint()[1], 0.0);
-        assertEquals(12.0, solution.getValue(), 0.0);
+        Assert.assertEquals(12.0, solution.getValue(), 0.0);
     }
 
     @Test(expected = NoFeasibleSolutionException.class)
@@ -203,12 +215,12 @@ public class SimplexSolverTest {
 
         SimplexSolver solver = new SimplexSolver();
         RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MAXIMIZE, true);
-        assertEquals(2902.92783505155, solution.getPoint()[0], .0000001);
+        Assert.assertEquals(2902.92783505155, solution.getPoint()[0], .0000001);
-        assertEquals(480.419243986254, solution.getPoint()[1], .0000001);
+        Assert.assertEquals(480.419243986254, solution.getPoint()[1], .0000001);
-        assertEquals(0.0, solution.getPoint()[2], .0000001);
+        Assert.assertEquals(0.0, solution.getPoint()[2], .0000001);
-        assertEquals(0.0, solution.getPoint()[3], .0000001);
+        Assert.assertEquals(0.0, solution.getPoint()[3], .0000001);
-        assertEquals(0.0, solution.getPoint()[4], .0000001);
+        Assert.assertEquals(0.0, solution.getPoint()[4], .0000001);
-        assertEquals(1438556.7491409, solution.getValue(), .0000001);
+        Assert.assertEquals(1438556.7491409, solution.getValue(), .0000001);
     }
 
     @Test
@@ -222,10 +234,10 @@ public class SimplexSolverTest {
 
       SimplexSolver solver = new SimplexSolver();
       RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MAXIMIZE, false);
-      assertEquals(1.0, solution.getPoint()[0], 0.0);
+      Assert.assertEquals(1.0, solution.getPoint()[0], 0.0);
-      assertEquals(1.0, solution.getPoint()[1], 0.0);
+      Assert.assertEquals(1.0, solution.getPoint()[1], 0.0);
-      assertEquals(0.0, solution.getPoint()[2], 0.0);
+      Assert.assertEquals(0.0, solution.getPoint()[2], 0.0);
-      assertEquals(15.0, solution.getValue(), 0.0);
+      Assert.assertEquals(15.0, solution.getValue(), 0.0);
   }
 
     @Test
@@ -236,7 +248,7 @@ public class SimplexSolverTest {
 
         SimplexSolver solver = new SimplexSolver();
         RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MAXIMIZE, true);
-        assertEquals(0, solution.getValue(), .0000001);
+        Assert.assertEquals(0, solution.getValue(), .0000001);
     }
 
     @Test
@@ -363,7 +375,7 @@ public class SimplexSolverTest {
 
         SimplexSolver solver = new SimplexSolver();
         RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MINIMIZE, true);
-        assertEquals(7518.0, solution.getValue(), .0000001);
+        Assert.assertEquals(7518.0, solution.getValue(), .0000001);
     }
 
     /**
@@ -385,13 +397,13 @@ public class SimplexSolverTest {
         String[] equationParts = s.split("[>|<]?=");
         double rhs = Double.parseDouble(equationParts[1].trim());
 
-        RealVector lhs = new ArrayRealVector(numCoefficients);
+        double[] lhs = new double[numCoefficients];
         String left = equationParts[0].replaceAll(" ?x", "");
         String[] coefficients = left.split(" ");
         for (String coefficient : coefficients) {
             double value = coefficient.charAt(0) == '-' ? -1 : 1;
             int index = Integer.parseInt(coefficient.replaceFirst("[+|-]", "").trim());
-            lhs.setEntry(index, value);
+            lhs[index] = value;
         }
         return new LinearConstraint(lhs, relationship, rhs);
     }

src/test/java/org/apache/commons/math/optimization/linear/SimplexTableauTest.java

@@ -22,11 +22,12 @@ import java.util.Collection;
 
 import org.apache.commons.math.TestUtils;
 import org.apache.commons.math.optimization.GoalType;
+import org.junit.Assert;
+import org.junit.Test;
 
-import junit.framework.TestCase;
+public class SimplexTableauTest {
-
-public class SimplexTableauTest extends TestCase {
 
+    @Test
     public void testInitialization() {
         LinearObjectiveFunction f = createFunction();
         Collection<LinearConstraint> constraints = createConstraints();
@@ -42,6 +43,7 @@ public class SimplexTableauTest extends TestCase {
         assertMatrixEquals(expectedInitialTableau, tableau.getData());
     }
 
+    @Test
     public void testdiscardArtificialVariables() {
         LinearObjectiveFunction f = createFunction();
         Collection<LinearConstraint> constraints = createConstraints();
@@ -57,6 +59,7 @@ public class SimplexTableauTest extends TestCase {
         assertMatrixEquals(expectedTableau, tableau.getData());
     }
 
+    @Test
     public void testTableauWithNoArtificialVars() {
         LinearObjectiveFunction f = new LinearObjectiveFunction(new double[] {15, 10}, 0);
         Collection<LinearConstraint> constraints = new ArrayList<LinearConstraint>();
@@ -74,13 +77,15 @@ public class SimplexTableauTest extends TestCase {
         assertMatrixEquals(initialTableau, tableau.getData());
     }
 
+    @Test
     public void testSerial() {
         LinearObjectiveFunction f = createFunction();
         Collection<LinearConstraint> constraints = createConstraints();
         SimplexTableau tableau =
             new SimplexTableau(f, constraints, GoalType.MAXIMIZE, false, 1.0e-6);
-        assertEquals(tableau, TestUtils.serializeAndRecover(tableau));
+        Assert.assertEquals(tableau, TestUtils.serializeAndRecover(tableau));
     }
+
     private LinearObjectiveFunction createFunction() {
         return new LinearObjectiveFunction(new double[] {15, 10}, 0);
     }
@@ -94,11 +99,11 @@ public class SimplexTableauTest extends TestCase {
     }
 
     private void assertMatrixEquals(double[][] expected, double[][] result) {
-        assertEquals("Wrong number of rows.", expected.length, result.length);
+        Assert.assertEquals("Wrong number of rows.", expected.length, result.length);
         for (int i = 0; i < expected.length; i++) {
-            assertEquals("Wrong number of columns.", expected[i].length, result[i].length);
+            Assert.assertEquals("Wrong number of columns.", expected[i].length, result[i].length);
             for (int j = 0; j < expected[i].length; j++) {
-                assertEquals("Wrong value at position [" + i + "," + j + "]", expected[i][j], result[i][j]);
+                Assert.assertEquals("Wrong value at position [" + i + "," + j + "]", expected[i][j], result[i][j], 1.0e-15);
             }
         }
     }