001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017
018package org.apache.commons.math3.optimization.linear;
019
020import java.util.ArrayList;
021import java.util.List;
022
023import org.apache.commons.math3.exception.MaxCountExceededException;
024import org.apache.commons.math3.optimization.PointValuePair;
025import org.apache.commons.math3.util.Precision;
026
027
028/**
029 * Solves a linear problem using the Two-Phase Simplex Method.
030 *
031 * @deprecated As of 3.1 (to be removed in 4.0).
032 * @since 2.0
033 */
034@Deprecated
035public class SimplexSolver extends AbstractLinearOptimizer {
036
037    /** Default amount of error to accept for algorithm convergence. */
038    private static final double DEFAULT_EPSILON = 1.0e-6;
039
040    /** Default amount of error to accept in floating point comparisons (as ulps). */
041    private static final int DEFAULT_ULPS = 10;
042
043    /** Amount of error to accept for algorithm convergence. */
044    private final double epsilon;
045
046    /** Amount of error to accept in floating point comparisons (as ulps). */
047    private final int maxUlps;
048
049    /**
050     * Build a simplex solver with default settings.
051     */
052    public SimplexSolver() {
053        this(DEFAULT_EPSILON, DEFAULT_ULPS);
054    }
055
056    /**
057     * Build a simplex solver with a specified accepted amount of error
058     * @param epsilon the amount of error to accept for algorithm convergence
059     * @param maxUlps amount of error to accept in floating point comparisons
060     */
061    public SimplexSolver(final double epsilon, final int maxUlps) {
062        this.epsilon = epsilon;
063        this.maxUlps = maxUlps;
064    }
065
066    /**
067     * Returns the column with the most negative coefficient in the objective function row.
068     * @param tableau simple tableau for the problem
069     * @return column with the most negative coefficient
070     */
071    private Integer getPivotColumn(SimplexTableau tableau) {
072        double minValue = 0;
073        Integer minPos = null;
074        for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getWidth() - 1; i++) {
075            final double entry = tableau.getEntry(0, i);
076            // check if the entry is strictly smaller than the current minimum
077            // do not use a ulp/epsilon check
078            if (entry < minValue) {
079                minValue = entry;
080                minPos = i;
081            }
082        }
083        return minPos;
084    }
085
086    /**
087     * Returns the row with the minimum ratio as given by the minimum ratio test (MRT).
088     * @param tableau simple tableau for the problem
089     * @param col the column to test the ratio of.  See {@link #getPivotColumn(SimplexTableau)}
090     * @return row with the minimum ratio
091     */
092    private Integer getPivotRow(SimplexTableau tableau, final int col) {
093        // create a list of all the rows that tie for the lowest score in the minimum ratio test
094        List<Integer> minRatioPositions = new ArrayList<Integer>();
095        double minRatio = Double.MAX_VALUE;
096        for (int i = tableau.getNumObjectiveFunctions(); i < tableau.getHeight(); i++) {
097            final double rhs = tableau.getEntry(i, tableau.getWidth() - 1);
098            final double entry = tableau.getEntry(i, col);
099
100            if (Precision.compareTo(entry, 0d, maxUlps) > 0) {
101                final double ratio = rhs / entry;
102                // check if the entry is strictly equal to the current min ratio
103                // do not use a ulp/epsilon check
104                final int cmp = Double.compare(ratio, minRatio);
105                if (cmp == 0) {
106                    minRatioPositions.add(i);
107                } else if (cmp < 0) {
108                    minRatio = ratio;
109                    minRatioPositions = new ArrayList<Integer>();
110                    minRatioPositions.add(i);
111                }
112            }
113        }
114
115        if (minRatioPositions.size() == 0) {
116            return null;
117        } else if (minRatioPositions.size() > 1) {
118            // there's a degeneracy as indicated by a tie in the minimum ratio test
119
120            // 1. check if there's an artificial variable that can be forced out of the basis
121            if (tableau.getNumArtificialVariables() > 0) {
122                for (Integer row : minRatioPositions) {
123                    for (int i = 0; i < tableau.getNumArtificialVariables(); i++) {
124                        int column = i + tableau.getArtificialVariableOffset();
125                        final double entry = tableau.getEntry(row, column);
126                        if (Precision.equals(entry, 1d, maxUlps) && row.equals(tableau.getBasicRow(column))) {
127                            return row;
128                        }
129                    }
130                }
131            }
132
133            // 2. apply Bland's rule to prevent cycling:
134            //    take the row for which the corresponding basic variable has the smallest index
135            //
136            // see http://www.stanford.edu/class/msande310/blandrule.pdf
137            // see http://en.wikipedia.org/wiki/Bland%27s_rule (not equivalent to the above paper)
138            //
139            // Additional heuristic: if we did not get a solution after half of maxIterations
140            //                       revert to the simple case of just returning the top-most row
141            // This heuristic is based on empirical data gathered while investigating MATH-828.
142            if (getIterations() < getMaxIterations() / 2) {
143                Integer minRow = null;
144                int minIndex = tableau.getWidth();
145                final int varStart = tableau.getNumObjectiveFunctions();
146                final int varEnd = tableau.getWidth() - 1;
147                for (Integer row : minRatioPositions) {
148                    for (int i = varStart; i < varEnd && !row.equals(minRow); i++) {
149                        final Integer basicRow = tableau.getBasicRow(i);
150                        if (basicRow != null && basicRow.equals(row) && i < minIndex) {
151                            minIndex = i;
152                            minRow = row;
153                        }
154                    }
155                }
156                return minRow;
157            }
158        }
159        return minRatioPositions.get(0);
160    }
161
162    /**
163     * Runs one iteration of the Simplex method on the given model.
164     * @param tableau simple tableau for the problem
165     * @throws MaxCountExceededException if the maximal iteration count has been exceeded
166     * @throws UnboundedSolutionException if the model is found not to have a bounded solution
167     */
168    protected void doIteration(final SimplexTableau tableau)
169        throws MaxCountExceededException, UnboundedSolutionException {
170
171        incrementIterationsCounter();
172
173        Integer pivotCol = getPivotColumn(tableau);
174        Integer pivotRow = getPivotRow(tableau, pivotCol);
175        if (pivotRow == null) {
176            throw new UnboundedSolutionException();
177        }
178
179        // set the pivot element to 1
180        double pivotVal = tableau.getEntry(pivotRow, pivotCol);
181        tableau.divideRow(pivotRow, pivotVal);
182
183        // set the rest of the pivot column to 0
184        for (int i = 0; i < tableau.getHeight(); i++) {
185            if (i != pivotRow) {
186                final double multiplier = tableau.getEntry(i, pivotCol);
187                tableau.subtractRow(i, pivotRow, multiplier);
188            }
189        }
190    }
191
192    /**
193     * Solves Phase 1 of the Simplex method.
194     * @param tableau simple tableau for the problem
195     * @throws MaxCountExceededException if the maximal iteration count has been exceeded
196     * @throws UnboundedSolutionException if the model is found not to have a bounded solution
197     * @throws NoFeasibleSolutionException if there is no feasible solution
198     */
199    protected void solvePhase1(final SimplexTableau tableau)
200        throws MaxCountExceededException, UnboundedSolutionException, NoFeasibleSolutionException {
201
202        // make sure we're in Phase 1
203        if (tableau.getNumArtificialVariables() == 0) {
204            return;
205        }
206
207        while (!tableau.isOptimal()) {
208            doIteration(tableau);
209        }
210
211        // if W is not zero then we have no feasible solution
212        if (!Precision.equals(tableau.getEntry(0, tableau.getRhsOffset()), 0d, epsilon)) {
213            throw new NoFeasibleSolutionException();
214        }
215    }
216
217    /** {@inheritDoc} */
218    @Override
219    public PointValuePair doOptimize()
220        throws MaxCountExceededException, UnboundedSolutionException, NoFeasibleSolutionException {
221        final SimplexTableau tableau =
222            new SimplexTableau(getFunction(),
223                               getConstraints(),
224                               getGoalType(),
225                               restrictToNonNegative(),
226                               epsilon,
227                               maxUlps);
228
229        solvePhase1(tableau);
230        tableau.dropPhase1Objective();
231
232        while (!tableau.isOptimal()) {
233            doIteration(tableau);
234        }
235        return tableau.getSolution();
236    }
237
238}