commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r1090656 - in /commons/proper/math/trunk: ./ src/main/java/org/apache/commons/math/optimization/linear/ src/site/xdoc/ src/test/java/org/apache/commons/math/optimization/linear/
Date Sat, 09 Apr 2011 19:20:48 GMT
Author: luc
Date: Sat Apr  9 19:20:47 2011
New Revision: 1090656

URL: http://svn.apache.org/viewvc?rev=1090656&view=rev
Log:
Fixed two errors in simplex solver when entries are close together or
when variables are not restricted to non-negative.

Jira: MATH-434

Modified:
    commons/proper/math/trunk/pom.xml
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/linear/SimplexSolver.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/linear/SimplexTableau.java
    commons/proper/math/trunk/src/site/xdoc/changes.xml
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/optimization/linear/SimplexSolverTest.java

Modified: commons/proper/math/trunk/pom.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/pom.xml?rev=1090656&r1=1090655&r2=1090656&view=diff
==============================================================================
--- commons/proper/math/trunk/pom.xml (original)
+++ commons/proper/math/trunk/pom.xml Sat Apr  9 19:20:47 2011
@@ -187,6 +187,9 @@
       <name>J. Lewis Muir</name>
     </contributor>
     <contributor>
+      <name>Thomas Neidhart</name>
+    </contributor>
+    <contributor>
       <name>Fredrik Norin</name>
     </contributor>
     <contributor>

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/linear/SimplexSolver.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/linear/SimplexSolver.java?rev=1090656&r1=1090655&r2=1090656&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/linear/SimplexSolver.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/linear/SimplexSolver.java
Sat Apr  9 19:20:47 2011
@@ -22,6 +22,7 @@ import java.util.List;
 
 import org.apache.commons.math.optimization.OptimizationException;
 import org.apache.commons.math.optimization.RealPointValuePair;
+import org.apache.commons.math.util.FastMath;
 import org.apache.commons.math.util.MathUtils;
 
 
@@ -31,26 +32,34 @@ import org.apache.commons.math.util.Math
  * @since 2.0
  */
 public class SimplexSolver extends AbstractLinearOptimizer {
-
-    /** Default amount of error to accept in floating point comparisons. */
+    
+    /** Default amount of error to accept for algorithm convergence. */
     private static final double DEFAULT_EPSILON = 1.0e-6;
-
-    /** Amount of error to accept in floating point comparisons. */
+     
+    /** Amount of error to accept for algorithm convergence. */
     protected final double epsilon;
 
+    /** Default amount of error to accept in floating point comparisons (as ulps). */
+    private static final int DEFAULT_ULPS = 10;
+
+    /** Amount of error to accept in floating point comparisons (as ulps). */
+    protected final int maxUlps;
+
     /**
      * Build a simplex solver with default settings.
      */
     public SimplexSolver() {
-        this(DEFAULT_EPSILON);
+        this(DEFAULT_EPSILON, DEFAULT_ULPS);
     }
 
     /**
      * Build a simplex solver with a specified accepted amount of error
-     * @param epsilon the amount of error to accept in floating point comparisons
+     * @param epsilon the amount of error to accept for algorithm convergence
+     * @param maxUlps amount of error to accept in floating point comparisons 
      */
-    public SimplexSolver(final double epsilon) {
+    public SimplexSolver(final double epsilon, final int maxUlps) {
         this.epsilon = epsilon;
+        this.maxUlps = maxUlps;
     }
 
     /**
@@ -62,8 +71,9 @@ public class SimplexSolver extends Abstr
         double minValue = 0;
         Integer minPos = null;
         for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getWidth() - 1; i++)
{
-            if (MathUtils.compareTo(tableau.getEntry(0, i), minValue, epsilon) < 0) {
-                minValue = tableau.getEntry(0, i);
+            final double entry = tableau.getEntry(0, i);
+            if (MathUtils.compareTo(entry, minValue, getEpsilon(entry)) < 0) {
+                minValue = entry;
                 minPos = i;
             }
         }
@@ -83,11 +93,13 @@ public class SimplexSolver extends Abstr
         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) {
+            
+            if (MathUtils.compareTo(entry, 0d, getEpsilon(entry)) > 0) {
                 final double ratio = rhs / entry;
-                if (MathUtils.equals(ratio, minRatio, epsilon)) {
+                final int cmp = MathUtils.compareTo(ratio, minRatio, getEpsilon(ratio));
+                if (cmp == 0) {
                     minRatioPositions.add(i);
-                } else if (ratio < minRatio) {
+                } else if (cmp < 0) {
                     minRatio = ratio;
                     minRatioPositions = new ArrayList<Integer>();
                     minRatioPositions.add(i);
@@ -103,7 +115,8 @@ public class SimplexSolver extends Abstr
           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) &&
+              final double entry = tableau.getEntry(row, column);
+              if (MathUtils.equals(entry, 1d, getEpsilon(entry)) &&
                   row.equals(tableau.getBasicRow(column))) {
                 return row;
               }
@@ -162,7 +175,7 @@ public class SimplexSolver extends Abstr
         }
 
         // if W is not zero then we have no feasible solution
-        if (!MathUtils.equals(tableau.getEntry(0, tableau.getRhsOffset()), 0, epsilon)) {
+        if (!MathUtils.equals(tableau.getEntry(0, tableau.getRhsOffset()), 0d, epsilon))
{
             throw new NoFeasibleSolutionException();
         }
     }
@@ -171,7 +184,8 @@ public class SimplexSolver extends Abstr
     @Override
     public RealPointValuePair doOptimize() throws OptimizationException {
         final SimplexTableau tableau =
-            new SimplexTableau(function, linearConstraints, goal, nonNegative, epsilon);
+            new SimplexTableau(function, linearConstraints, goal, nonNegative, 
+                               epsilon, maxUlps);
 
         solvePhase1(tableau);
         tableau.dropPhase1Objective();
@@ -182,4 +196,12 @@ public class SimplexSolver extends Abstr
         return tableau.getSolution();
     }
 
+    /**
+     * Get an epsilon that is adjusted to the magnitude of the given value.
+     * @param value the value for which to get the epsilon
+     * @return magnitude-adjusted epsilon using {@link FastMath.ulp}
+     */
+    private double getEpsilon(double value) {
+        return FastMath.ulp(value) * (double) maxUlps;
+    }
 }

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/linear/SimplexTableau.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/linear/SimplexTableau.java?rev=1090656&r1=1090655&r2=1090656&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/linear/SimplexTableau.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/optimization/linear/SimplexTableau.java
Sat Apr  9 19:20:47 2011
@@ -33,6 +33,7 @@ import org.apache.commons.math.linear.Re
 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.FastMath;
 import org.apache.commons.math.util.MathUtils;
 
 /**
@@ -65,6 +66,9 @@ class SimplexTableau implements Serializ
     /** Column label for negative vars. */
     private static final String NEGATIVE_VAR_COLUMN_LABEL = "x-";
 
+    /** Default amount of error to accept in floating point comparisons (as ulps). */
+    private static final int DEFAULT_ULPS = 10;
+
     /** Serializable version identifier. */
     private static final long serialVersionUID = -1369660067587938365L;
 
@@ -92,9 +96,12 @@ class SimplexTableau implements Serializ
     /** Number of artificial variables. */
     private int numArtificialVariables;
 
-    /** Amount of error to accept in floating point comparisons. */
+    /** Amount of error to accept when checking for optimality. */
     private final double epsilon;
 
+    /** Amount of error to accept in floating point comparisons. */
+    private final int maxUlps;
+
     /**
      * Build a tableau for a linear problem.
      * @param f linear objective function
@@ -102,16 +109,35 @@ class SimplexTableau implements Serializ
      * @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
+     * @param epsilon amount of error to accept when checking for optimality
      */
     SimplexTableau(final LinearObjectiveFunction f,
                    final Collection<LinearConstraint> constraints,
                    final GoalType goalType, final boolean restrictToNonNegative,
                    final double epsilon) {
+        this(f, constraints, goalType, restrictToNonNegative, epsilon, DEFAULT_ULPS);
+    }
+    
+    /**
+     * 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 when checking for optimality
+     * @param maxUlps 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,
+                   final int maxUlps) {
         this.f                      = f;
         this.constraints            = normalizeConstraints(constraints);
         this.restrictToNonNegative  = restrictToNonNegative;
         this.epsilon                = epsilon;
+        this.maxUlps                = maxUlps;
         this.numDecisionVariables   = f.getCoefficients().getDimension() +
                                       (restrictToNonNegative ? 0 : 1);
         this.numSlackVariables      = getConstraintTypeCounts(Relationship.LEQ) +
@@ -172,7 +198,7 @@ class SimplexTableau implements Serializ
 
         if (!restrictToNonNegative) {
             matrix.setEntry(zIndex, getSlackVariableOffset() - 1,
-                getInvertedCoeffiecientSum(objectiveCoefficients));
+                getInvertedCoefficientSum(objectiveCoefficients));
         }
 
         // initialize the constraint rows
@@ -188,7 +214,7 @@ class SimplexTableau implements Serializ
             // x-
             if (!restrictToNonNegative) {
                 matrix.setEntry(row, getSlackVariableOffset() - 1,
-                    getInvertedCoeffiecientSum(constraint.getCoefficients()));
+                    getInvertedCoefficientSum(constraint.getCoefficients()));
             }
 
             // RHS
@@ -269,7 +295,7 @@ class SimplexTableau implements Serializ
      * @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) {
+    protected static double getInvertedCoefficientSum(final RealVector coefficients) {
         double sum = 0;
         for (double coefficient : coefficients.getData()) {
             sum -= coefficient;
@@ -285,9 +311,10 @@ class SimplexTableau implements Serializ
     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))
{
+            final double entry = getEntry(i, col);
+            if (MathUtils.equals(entry, 1d, getEpsilon(entry)) && (row == null))
{
                 row = i;
-            } else if (!MathUtils.equals(getEntry(i, col), 0.0, epsilon)) {
+            } else if (!MathUtils.equals(entry, 0d, getEpsilon(entry))) {
                 return null;
             }
         }
@@ -308,9 +335,10 @@ class SimplexTableau implements Serializ
 
         // 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);
-          }
+            final double entry = tableau.getEntry(0, i);
+            if (MathUtils.compareTo(entry, 0d, getEpsilon(entry)) > 0) {
+                columnsToDrop.add(i);
+            }
         }
 
         // non-basic artificial variables
@@ -353,7 +381,8 @@ class SimplexTableau implements Serializ
      */
     boolean isOptimal() {
         for (int i = getNumObjectiveFunctions(); i < getWidth() - 1; i++) {
-            if (MathUtils.compareTo(tableau.getEntry(0, i), 0, epsilon) < 0) {
+            final double entry = tableau.getEntry(0, i);
+            if (MathUtils.compareTo(entry, 0d, epsilon) < 0) {
                 return false;
             }
         }
@@ -382,7 +411,7 @@ class SimplexTableau implements Serializ
           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;
+              coefficients[i] = 0 - (restrictToNonNegative ? 0 : mostNegative);
           } else {
               basicRows.add(basicRow);
               coefficients[i] =
@@ -545,6 +574,7 @@ class SimplexTableau implements Serializ
                  (numSlackVariables      == rhs.numSlackVariables) &&
                  (numArtificialVariables == rhs.numArtificialVariables) &&
                  (epsilon                == rhs.epsilon) &&
+                 (maxUlps                == rhs.maxUlps) &&
                  f.equals(rhs.f) &&
                  constraints.equals(rhs.constraints) &&
                  tableau.equals(rhs.tableau);
@@ -560,6 +590,7 @@ class SimplexTableau implements Serializ
                numSlackVariables ^
                numArtificialVariables ^
                Double.valueOf(epsilon).hashCode() ^
+               maxUlps ^
                f.hashCode() ^
                constraints.hashCode() ^
                tableau.hashCode();
@@ -585,5 +616,13 @@ class SimplexTableau implements Serializ
         ois.defaultReadObject();
         MatrixUtils.deserializeRealMatrix(this, "tableau", ois);
     }
-
+    
+    /**
+     * Get an epsilon that is adjusted to the magnitude of the given value.
+     * @param value the value for which to get the epsilon
+     * @return magnitude-adjusted epsilon using {@link FastMath.ulp}
+     */
+    private double getEpsilon(double value) {
+        return FastMath.ulp(value) * (double) maxUlps;
+    }    
 }

Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=1090656&r1=1090655&r2=1090656&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Sat Apr  9 19:20:47 2011
@@ -52,6 +52,10 @@ The <action> type attribute can be add,u
     If the output is not quite correct, check for invisible trailing spaces!
      -->
     <release version="3.0" date="TBD" description="TBD">
+      <action dev="luc" type="fix" issue="MATH-434" due-to="Thomas Neidhart">
+        Fixed two errors in simplex solver when entries are close together or
+        when variables are not restricted to non-negative.
+      </action>
       <action dev="luc" type="fix" issue="MATH-547" due-to="Thomas Neidhart">
         Improved robustness of k-means++ algorithm, by tracking changes in points assignments
         to clusters.

Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/optimization/linear/SimplexSolverTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/optimization/linear/SimplexSolverTest.java?rev=1090656&r1=1090655&r2=1090656&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/optimization/linear/SimplexSolverTest.java
(original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/optimization/linear/SimplexSolverTest.java
Sat Apr  9 19:20:47 2011
@@ -25,11 +25,88 @@ import java.util.Collection;
 import org.apache.commons.math.optimization.GoalType;
 import org.apache.commons.math.optimization.OptimizationException;
 import org.apache.commons.math.optimization.RealPointValuePair;
+import org.apache.commons.math.util.MathUtils;
 import org.junit.Test;
 
 public class SimplexSolverTest {
 
     @Test
+    public void test434NegativeVariable() throws OptimizationException
+    {
+        LinearObjectiveFunction f = new LinearObjectiveFunction(new double[] {0.0, 0.0, 1.0},
0.0d);
+        ArrayList<LinearConstraint> constraints = new ArrayList<LinearConstraint>();
+        constraints.add(new LinearConstraint(new double[] {1, 1, 0}, Relationship.EQ, 5));
+        constraints.add(new LinearConstraint(new double[] {0, 0, 1}, Relationship.GEQ, -10));
+
+        double epsilon = 1e-6;
+        SimplexSolver solver = new SimplexSolver();
+        RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MINIMIZE,
false);
+
+        Assert.assertEquals(5.0, solution.getPoint()[0] + solution.getPoint()[1], epsilon);
+        Assert.assertEquals(-10.0, solution.getPoint()[2], epsilon);
+        Assert.assertEquals(-10.0, solution.getValue(), epsilon);
+
+    }
+
+    @Test(expected = NoFeasibleSolutionException.class)
+    public void test434UnfeasibleSolution() throws OptimizationException
+    {
+        double epsilon = 1e-6;
+
+        LinearObjectiveFunction f = new LinearObjectiveFunction(new double[] {1.0, 0.0},
0.0);
+        ArrayList<LinearConstraint> constraints = new ArrayList<LinearConstraint>();
+        constraints.add(new LinearConstraint(new double[] {epsilon/2, 0.5}, Relationship.EQ,
0));
+        constraints.add(new LinearConstraint(new double[] {1e-3, 0.1}, Relationship.EQ, 10));
+
+        SimplexSolver solver = new SimplexSolver();
+        // allowing only non-negative values, no feasible solution shall be found
+        solver.optimize(f, constraints, GoalType.MINIMIZE, true);
+    }
+
+    @Test
+    public void test434PivotRowSelection() throws OptimizationException
+    {
+        LinearObjectiveFunction f = new LinearObjectiveFunction(new double[] {1.0}, 0.0);
+
+        double epsilon = 1e-6;
+        ArrayList<LinearConstraint> constraints = new ArrayList<LinearConstraint>();
+        constraints.add(new LinearConstraint(new double[] {200}, Relationship.GEQ, 1));
+        constraints.add(new LinearConstraint(new double[] {100}, Relationship.GEQ, 0.499900001));
+
+        SimplexSolver solver = new SimplexSolver();
+        RealPointValuePair solution = solver.optimize(f, constraints, GoalType.MINIMIZE,
false);
+        
+        Assert.assertTrue(MathUtils.compareTo(solution.getPoint()[0] * 200.d, 1.d, epsilon)
>= 0);
+        Assert.assertEquals(0.0050, solution.getValue(), epsilon);
+    }
+
+    @Test
+    public void test434PivotRowSelection2() throws OptimizationException
+    {
+        LinearObjectiveFunction f = new LinearObjectiveFunction(new double[] {0.0d, 1.0d,
1.0d, 0.0d, 0.0d, 0.0d, 0.0d}, 0.0d);
+
+        ArrayList<LinearConstraint> constraints = new ArrayList<LinearConstraint>();
+        constraints.add(new LinearConstraint(new double[] {1.0d, -0.1d, 0.0d, 0.0d, 0.0d,
0.0d, 0.0d}, Relationship.EQ, -0.1d));
+        constraints.add(new LinearConstraint(new double[] {1.0d, 0.0d, 0.0d, 0.0d, 0.0d,
0.0d, 0.0d}, Relationship.GEQ, -1e-18d));
+        constraints.add(new LinearConstraint(new double[] {0.0d, 1.0d, 0.0d, 0.0d, 0.0d,
0.0d, 0.0d}, Relationship.GEQ, 0.0d));
+        constraints.add(new LinearConstraint(new double[] {0.0d, 0.0d, 0.0d, 1.0d, 0.0d,
-0.0128588d, 1e-5d}, Relationship.EQ, 0.0d));
+        constraints.add(new LinearConstraint(new double[] {0.0d, 0.0d, 0.0d, 0.0d, 1.0d,
1e-5d, -0.0128586d}, Relationship.EQ, 1e-10d));
+        constraints.add(new LinearConstraint(new double[] {0.0d, 0.0d, 1.0d, -1.0d, 0.0d,
0.0d, 0.0d}, Relationship.GEQ, 0.0d));
+        constraints.add(new LinearConstraint(new double[] {0.0d, 0.0d, 1.0d, 1.0d, 0.0d,
0.0d, 0.0d}, Relationship.GEQ, 0.0d));
+        constraints.add(new LinearConstraint(new double[] {0.0d, 0.0d, 1.0d, 0.0d, -1.0d,
0.0d, 0.0d}, Relationship.GEQ, 0.0d));
+        constraints.add(new LinearConstraint(new double[] {0.0d, 0.0d, 1.0d, 0.0d, 1.0d,
0.0d, 0.0d}, Relationship.GEQ, 0.0d));
+
+        double epsilon = 1e-7;
+        SimplexSolver simplex = new SimplexSolver();
+        RealPointValuePair solution = simplex.optimize(f, constraints, GoalType.MINIMIZE,
false);
+        
+        Assert.assertTrue(MathUtils.compareTo(solution.getPoint()[0], -1e-18d, epsilon) >=
0);
+        Assert.assertEquals(1.0d, solution.getPoint()[1], epsilon);        
+        Assert.assertEquals(0.0d, solution.getPoint()[2], epsilon);
+        Assert.assertEquals(1.0d, solution.getValue(), epsilon);
+    }
+    
+    @Test
     public void testMath272() throws OptimizationException {
         LinearObjectiveFunction f = new LinearObjectiveFunction(new double[] { 2, 2, 1 },
0);
         Collection<LinearConstraint> constraints = new ArrayList<LinearConstraint>();



Mime
View raw message