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>();
