SimplexTableau.java

  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. package org.apache.commons.math4.legacy.optim.linear;

  18. import java.util.ArrayList;
  19. import java.util.Arrays;
  20. import java.util.Collection;
  21. import java.util.HashSet;
  22. import java.util.List;
  23. import java.util.Set;
  24. import java.util.TreeSet;

  25. import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
  26. import org.apache.commons.math4.legacy.linear.Array2DRowRealMatrix;
  27. import org.apache.commons.math4.legacy.linear.RealVector;
  28. import org.apache.commons.math4.legacy.optim.PointValuePair;
  29. import org.apache.commons.math4.legacy.optim.nonlinear.scalar.GoalType;
  30. import org.apache.commons.numbers.core.Precision;

  31. /**
  32.  * A tableau for use in the Simplex method.
  33.  *
  34.  * <p>
  35.  * Example:
  36.  * <pre>
  37.  *   W |  Z |  x1 |  x2 |  x- | s1 |  s2 |  a1 |  RHS
  38.  * ---------------------------------------------------
  39.  *  -1    0    0     0     0     0     0     1     0   &lt;= phase 1 objective
  40.  *   0    1   -15   -10    0     0     0     0     0   &lt;= phase 2 objective
  41.  *   0    0    1     0     0     1     0     0     2   &lt;= constraint 1
  42.  *   0    0    0     1     0     0     1     0     3   &lt;= constraint 2
  43.  *   0    0    1     1     0     0     0     1     4   &lt;= constraint 3
  44.  * </pre>
  45.  * W: Phase 1 objective function<br>
  46.  * Z: Phase 2 objective function<br>
  47.  * x1 &amp; x2: Decision variables<br>
  48.  * x-: Extra decision variable to allow for negative values<br>
  49.  * s1 &amp; s2: Slack/Surplus variables<br>
  50.  * a1: Artificial variable<br>
  51.  * RHS: Right hand side<br>
  52.  *
  53.  * Note on usage and safety:
  54.  * The class is package private. It is not meant for public usage.
  55.  * The core data structure, the tableau field, is mutated internally and
  56.  * even reallocated when necessary.
  57.  * Proper usage of this class is demonstrated in SimplexSolver,
  58.  * where the class is only ever constructed in a method (never a field
  59.  * of an object), and its lifetime, is therefore bound to a single thread (the
  60.  * thread that's invoking the method).
  61.  *
  62.  * @since 2.0
  63.  */
  64. class SimplexTableau {

  65.     /** Column label for negative vars. */
  66.     private static final String NEGATIVE_VAR_COLUMN_LABEL = "x-";
  67.     /** bit mask for IEEE double exponent. */
  68.     private static final long EXPN = 0x7ff0000000000000L;
  69.     /** bit mask for IEEE double mantissa and sign. */
  70.     private static final long FRAC = 0x800fffffffffffffL;
  71.     /** max IEEE exponent is 2047. */
  72.     private static final int MAX_IEEE_EXP = 2047;
  73.     /** min IEEE exponent is 0. */
  74.     private static final int MIN_IEEE_EXP = 0;
  75.     /** IEEE exponent is kept in an offset form, 1023 is zero. */
  76.     private static final int OFFSET_IEEE_EXP = 1023;
  77.     /** double exponent shift per IEEE standard. */
  78.     private static final int IEEE_EXPONENT_SHIFT = 52;

  79.     /** Linear objective function. */
  80.     private final LinearObjectiveFunction f;

  81.     /** Linear constraints. */
  82.     private final List<LinearConstraint> constraints;

  83.     /** Whether to restrict the variables to non-negative values. */
  84.     private final boolean restrictToNonNegative;

  85.     /** The variables each column represents. */
  86.     private final List<String> columnLabels = new ArrayList<>();

  87.     /** Simple tableau. */
  88.     private Array2DRowRealMatrix tableau;

  89.     /** Number of decision variables. */
  90.     private final int numDecisionVariables;

  91.     /** Number of slack variables. */
  92.     private final int numSlackVariables;

  93.     /** Number of artificial variables. */
  94.     private int numArtificialVariables;

  95.     /** Amount of error to accept when checking for optimality. */
  96.     private final double epsilon;

  97.     /** Amount of error to accept in floating point comparisons. */
  98.     private final int maxUlps;

  99.     /** Maps basic variables to row they are basic in. */
  100.     private int[] basicVariables;

  101.     /** Maps rows to their corresponding basic variables. */
  102.     private int[] basicRows;

  103.     /** changes in floating point exponent to scale the input. */
  104.     private int[] variableExpChange;

  105.     /**
  106.      * Builds a tableau for a linear problem.
  107.      *
  108.      * @param f Linear objective function.
  109.      * @param constraints Linear constraints.
  110.      * @param goalType Optimization goal: either {@link GoalType#MAXIMIZE}
  111.      * or {@link GoalType#MINIMIZE}.
  112.      * @param restrictToNonNegative Whether to restrict the variables to non-negative values.
  113.      * @param epsilon Amount of error to accept when checking for optimality.
  114.      * @throws DimensionMismatchException if the dimension of the constraints does not match the
  115.      *   dimension of the objective function
  116.      */
  117.     SimplexTableau(final LinearObjectiveFunction f,
  118.                    final Collection<LinearConstraint> constraints,
  119.                    final GoalType goalType,
  120.                    final boolean restrictToNonNegative,
  121.                    final double epsilon) {
  122.         this(f, constraints, goalType, restrictToNonNegative, epsilon, SimplexSolver.DEFAULT_ULPS);
  123.     }

  124.     /**
  125.      * Build a tableau for a linear problem.
  126.      * @param f linear objective function
  127.      * @param constraints linear constraints
  128.      * @param goalType type of optimization goal: either {@link GoalType#MAXIMIZE} or {@link GoalType#MINIMIZE}
  129.      * @param restrictToNonNegative whether to restrict the variables to non-negative values
  130.      * @param epsilon amount of error to accept when checking for optimality
  131.      * @param maxUlps amount of error to accept in floating point comparisons
  132.      * @throws DimensionMismatchException if the dimension of the constraints does not match the
  133.      *   dimension of the objective function
  134.      */
  135.     SimplexTableau(final LinearObjectiveFunction f,
  136.                    final Collection<LinearConstraint> constraints,
  137.                    final GoalType goalType,
  138.                    final boolean restrictToNonNegative,
  139.                    final double epsilon,
  140.                    final int maxUlps) throws DimensionMismatchException {
  141.         checkDimensions(f, constraints);
  142.         this.f                      = f;
  143.         this.constraints            = normalizeConstraints(constraints);
  144.         this.restrictToNonNegative  = restrictToNonNegative;
  145.         this.epsilon                = epsilon;
  146.         this.maxUlps                = maxUlps;
  147.         this.numDecisionVariables   = f.getCoefficients().getDimension() + (restrictToNonNegative ? 0 : 1);
  148.         this.numSlackVariables      = getConstraintTypeCounts(Relationship.LEQ) +
  149.                                       getConstraintTypeCounts(Relationship.GEQ);
  150.         this.numArtificialVariables = getConstraintTypeCounts(Relationship.EQ) +
  151.                                       getConstraintTypeCounts(Relationship.GEQ);
  152.         this.tableau = createTableau(goalType == GoalType.MAXIMIZE);
  153.         // initialize the basic variables for phase 1:
  154.         //   we know that only slack or artificial variables can be basic
  155.         initializeBasicVariables(getSlackVariableOffset());
  156.         initializeColumnLabels();
  157.     }

  158.     /**
  159.      * Checks that the dimensions of the objective function and the constraints match.
  160.      * @param objectiveFunction the objective function
  161.      * @param c the set of constraints
  162.      * @throws DimensionMismatchException if the constraint dimensions do not match with the
  163.      *   dimension of the objective function
  164.      */
  165.     private void checkDimensions(final LinearObjectiveFunction objectiveFunction,
  166.                                  final Collection<LinearConstraint> c) {
  167.         final int dimension = objectiveFunction.getCoefficients().getDimension();
  168.         for (final LinearConstraint constraint : c) {
  169.             final int constraintDimension = constraint.getCoefficients().getDimension();
  170.             if (constraintDimension != dimension) {
  171.                 throw new DimensionMismatchException(constraintDimension, dimension);
  172.             }
  173.         }
  174.     }
  175.     /**
  176.      * Initialize the labels for the columns.
  177.      */
  178.     protected void initializeColumnLabels() {
  179.       if (getNumObjectiveFunctions() == 2) {
  180.         columnLabels.add("W");
  181.       }
  182.       columnLabels.add("Z");
  183.       for (int i = 0; i < getOriginalNumDecisionVariables(); i++) {
  184.         columnLabels.add("x" + i);
  185.       }
  186.       if (!restrictToNonNegative) {
  187.         columnLabels.add(NEGATIVE_VAR_COLUMN_LABEL);
  188.       }
  189.       for (int i = 0; i < getNumSlackVariables(); i++) {
  190.         columnLabels.add("s" + i);
  191.       }
  192.       for (int i = 0; i < getNumArtificialVariables(); i++) {
  193.         columnLabels.add("a" + i);
  194.       }
  195.       columnLabels.add("RHS");
  196.     }

  197.     /**
  198.      * Create the tableau by itself.
  199.      * @param maximize if true, goal is to maximize the objective function
  200.      * @return created tableau
  201.      */
  202.     protected Array2DRowRealMatrix createTableau(final boolean maximize) {

  203.         // create a matrix of the correct size
  204.         int width = numDecisionVariables + numSlackVariables +
  205.         numArtificialVariables + getNumObjectiveFunctions() + 1; // + 1 is for RHS
  206.         int height = constraints.size() + getNumObjectiveFunctions();
  207.         Array2DRowRealMatrix matrix = new Array2DRowRealMatrix(height, width);

  208.         // initialize the objective function rows
  209.         if (getNumObjectiveFunctions() == 2) {
  210.             matrix.setEntry(0, 0, -1);
  211.         }

  212.         int zIndex = (getNumObjectiveFunctions() == 1) ? 0 : 1;
  213.         matrix.setEntry(zIndex, zIndex, maximize ? 1 : -1);

  214.         double[][] scaled = new double[constraints.size() + 1][];

  215.         RealVector objectiveCoefficients = maximize ? f.getCoefficients().mapMultiply(-1) : f.getCoefficients();
  216.         scaled[0] = objectiveCoefficients.toArray();
  217.         double[] scaledRhs = new double[constraints.size() + 1];
  218.         double value = maximize ? f.getConstantTerm() : -1 * f.getConstantTerm();
  219.         scaledRhs[0] = value;

  220.         for (int i = 0; i < constraints.size(); i++) {
  221.             LinearConstraint constraint = constraints.get(i);
  222.             scaled[i + 1] = constraint.getCoefficients().toArray();
  223.             scaledRhs[i + 1] = constraint.getValue();
  224.         }
  225.         variableExpChange = new int[scaled[0].length];

  226.         scale(scaled, scaledRhs);

  227.         copyArray(scaled[0], matrix.getDataRef()[zIndex]);
  228.         matrix.setEntry(zIndex, width - 1, scaledRhs[0]);

  229.         if (!restrictToNonNegative) {
  230.             matrix.setEntry(zIndex, getSlackVariableOffset() - 1,
  231.                             getInvertedCoefficientSum(scaled[0]));
  232.         }

  233.         // initialize the constraint rows
  234.         int slackVar = 0;
  235.         int artificialVar = 0;
  236.         for (int i = 0; i < constraints.size(); i++) {
  237.             final LinearConstraint constraint = constraints.get(i);
  238.             final int row = getNumObjectiveFunctions() + i;

  239.             // decision variable coefficients
  240.             copyArray(scaled[i + 1], matrix.getDataRef()[row]);

  241.             // x-
  242.             if (!restrictToNonNegative) {
  243.                 matrix.setEntry(row, getSlackVariableOffset() - 1,
  244.                                 getInvertedCoefficientSum(scaled[i + 1]));
  245.             }

  246.             // RHS
  247.             matrix.setEntry(row, width - 1, scaledRhs[i + 1]);

  248.             // slack variables
  249.             if (constraint.getRelationship() == Relationship.LEQ) {
  250.                 matrix.setEntry(row, getSlackVariableOffset() + slackVar++, 1);  // slack
  251.             } else if (constraint.getRelationship() == Relationship.GEQ) {
  252.                 matrix.setEntry(row, getSlackVariableOffset() + slackVar++, -1); // excess
  253.             }

  254.             // artificial variables
  255.             if (constraint.getRelationship() == Relationship.EQ ||
  256.                 constraint.getRelationship() == Relationship.GEQ) {
  257.                 matrix.setEntry(0, getArtificialVariableOffset() + artificialVar, 1);
  258.                 matrix.setEntry(row, getArtificialVariableOffset() + artificialVar++, 1);
  259.                 matrix.setRowVector(0, matrix.getRowVector(0).subtract(matrix.getRowVector(row)));
  260.             }
  261.         }

  262.         return matrix;
  263.     }

  264.     /** We scale the constants in the equations and objective, which means we try
  265.      * to get the IEEE double exponent as close to zero (1023) as possible, which makes the
  266.      * constants closer to 1.
  267.      * We use exponent shifts instead of division because that introduces no bit errors.
  268.      *
  269.      * @param scaled coefficients before scaling
  270.      * @param scaledRhs right hand side before scaling
  271.      */
  272.     private void scale(double[][] scaled, double[] scaledRhs) {
  273.         /*
  274.             first transform across:
  275.             c0 x0 + c1 x1 + ... + cn xn = vn ==> (2^expChange) * (c0 x0 + c1 x1 + ... + cn xn = vn)

  276.             expChange will be negative if the constants are larger than 1,
  277.             it'll be positive if the constants are less than 1.
  278.         */
  279.         for (int i = 0; i < scaled.length; i++) {
  280.             int minExp = MAX_IEEE_EXP + 1;
  281.             int maxExp = MIN_IEEE_EXP - 1;
  282.             for (double d: scaled[i]) {
  283.                 if (d != 0) {
  284.                     int e = exponent(d);
  285.                     if (e < minExp) {
  286.                         minExp = e;
  287.                     }
  288.                     if (e > maxExp) {
  289.                         maxExp = e;
  290.                     }
  291.                 }
  292.             }
  293.             if (scaledRhs[i] != 0) {
  294.                 final int e = exponent(scaledRhs[i]);
  295.                 if (e < minExp) {
  296.                     minExp = e;
  297.                 }
  298.                 if (e > maxExp) {
  299.                     maxExp = e;
  300.                 }
  301.             }
  302.             final int expChange = computeExpChange(minExp, maxExp);
  303.             if (expChange != 0) {
  304.                 scaledRhs[i] = updateExponent(scaledRhs[i], expChange);
  305.                 updateExponent(scaled[i], expChange);
  306.             }
  307.         }

  308.         /*
  309.             second, transform down the columns. this is like defining a new variable for that column
  310.             that is yi = xi * (2^expChange)
  311.             After solving for yi, we compute xi by shifting again. See getSolution()
  312.          */
  313.         for (int i = 0; i < variableExpChange.length; i++) {
  314.             int minExp = MAX_IEEE_EXP + 1;
  315.             int maxExp = MIN_IEEE_EXP - 1;

  316.             for (double[] coefficients : scaled) {
  317.                 final double d = coefficients[i];
  318.                 if (d != 0) {
  319.                     int e = exponent(d);
  320.                     if (e < minExp) {
  321.                         minExp = e;
  322.                     }
  323.                     if (e > maxExp) {
  324.                         maxExp = e;
  325.                     }
  326.                 }
  327.             }
  328.             final int expChange = computeExpChange(minExp, maxExp);
  329.             variableExpChange[i] = expChange;
  330.             if (expChange != 0) {
  331.                 for (double[] coefficients : scaled) {
  332.                      coefficients[i] = updateExponent(coefficients[i], expChange);
  333.                 }
  334.             }
  335.         }
  336.     }

  337.     /**
  338.      * Given the minimum and maximum value of the exponent of two {@code double}
  339.      * values, pick a change in exponent to bring those values closer to 1.
  340.      *
  341.      * @param minExp Smallest exponent.
  342.      * @param maxExp Largest exponent.
  343.      * @return the new exponent.
  344.      */
  345.     private int computeExpChange(int minExp, int maxExp) {
  346.         int expChange = 0;
  347.         if (minExp <= MAX_IEEE_EXP &&
  348.             minExp > OFFSET_IEEE_EXP) {
  349.             expChange = OFFSET_IEEE_EXP - minExp;
  350.         } else if (maxExp >= MIN_IEEE_EXP &&
  351.                    maxExp < OFFSET_IEEE_EXP) {
  352.             expChange = OFFSET_IEEE_EXP - maxExp;
  353.         }
  354.         return expChange;
  355.     }

  356.     /**
  357.      * Changes the exponent of every member of the array by the given amount.
  358.      *
  359.      * @param dar array of doubles to change
  360.      * @param exp exponent value to change
  361.      */
  362.     private static void updateExponent(double[] dar, int exp) {
  363.         for (int i = 0; i < dar.length; i++) {
  364.             dar[i] = updateExponent(dar[i], exp);
  365.         }
  366.     }

  367.     /**
  368.      * Extract the exponent of a {@code double}.
  369.      *
  370.      * @param d value to extract the exponent from
  371.      * @return the IEEE exponent in the EXPN bits, as an integer
  372.      */
  373.     private static int exponent(double d) {
  374.         final long bits = Double.doubleToLongBits(d);
  375.         return (int) ((bits & EXPN) >>> IEEE_EXPONENT_SHIFT);
  376.     }

  377.     /**
  378.      * Changes the exponent of a number by the given amount.
  379.      *
  380.      * @param d value to change
  381.      * @param exp exponent to add to the existing exponent (may be negative)
  382.      * @return a double with the same sign/mantissa bits as d, but exponent changed by exp
  383.      */
  384.     private static double updateExponent(double d, int exp) {
  385.         if (d == 0 ||
  386.             exp == 0) {
  387.             return d;
  388.         }
  389.         final long bits = Double.doubleToLongBits(d);
  390.         return Double.longBitsToDouble((bits & FRAC) | ((((bits & EXPN) >>> IEEE_EXPONENT_SHIFT) + exp) << IEEE_EXPONENT_SHIFT));
  391.     }

  392.     /**
  393.      * Get new versions of the constraints which have positive right hand sides.
  394.      * @param originalConstraints original (not normalized) constraints
  395.      * @return new versions of the constraints
  396.      */
  397.     public List<LinearConstraint> normalizeConstraints(Collection<LinearConstraint> originalConstraints) {
  398.         final List<LinearConstraint> normalized = new ArrayList<>(originalConstraints.size());
  399.         for (LinearConstraint constraint : originalConstraints) {
  400.             normalized.add(normalize(constraint));
  401.         }
  402.         return normalized;
  403.     }

  404.     /**
  405.      * Get a new equation equivalent to this one with a positive right hand side.
  406.      * @param constraint reference constraint
  407.      * @return new equation
  408.      */
  409.     private LinearConstraint normalize(final LinearConstraint constraint) {
  410.         if (constraint.getValue() < 0) {
  411.             return new LinearConstraint(constraint.getCoefficients().mapMultiply(-1),
  412.                                         constraint.getRelationship().oppositeRelationship(),
  413.                                         -1 * constraint.getValue());
  414.         }
  415.         return new LinearConstraint(constraint.getCoefficients(),
  416.                                     constraint.getRelationship(), constraint.getValue());
  417.     }

  418.     /**
  419.      * Get the number of objective functions in this tableau.
  420.      * @return 2 for Phase 1.  1 for Phase 2.
  421.      */
  422.     protected final int getNumObjectiveFunctions() {
  423.         return this.numArtificialVariables > 0 ? 2 : 1;
  424.     }

  425.     /**
  426.      * Get a count of constraints corresponding to a specified relationship.
  427.      * @param relationship relationship to count
  428.      * @return number of constraint with the specified relationship
  429.      */
  430.     private int getConstraintTypeCounts(final Relationship relationship) {
  431.         int count = 0;
  432.         for (final LinearConstraint constraint : constraints) {
  433.             if (constraint.getRelationship() == relationship) {
  434.                 ++count;
  435.             }
  436.         }
  437.         return count;
  438.     }

  439.     /**
  440.      * Get the -1 times the sum of all coefficients in the given array.
  441.      * @param coefficients coefficients to sum
  442.      * @return the -1 times the sum of all coefficients in the given array.
  443.      */
  444.     private static double getInvertedCoefficientSum(final double[] coefficients) {
  445.         double sum = 0;
  446.         for (double coefficient : coefficients) {
  447.             sum -= coefficient;
  448.         }
  449.         return sum;
  450.     }

  451.     /**
  452.      * Checks whether the given column is basic.
  453.      * @param col index of the column to check
  454.      * @return the row that the variable is basic in.  null if the column is not basic
  455.      */
  456.     protected Integer getBasicRow(final int col) {
  457.         final int row = basicVariables[col];
  458.         return row == -1 ? null : row;
  459.     }

  460.     /**
  461.      * Returns the variable that is basic in this row.
  462.      * @param row the index of the row to check
  463.      * @return the variable that is basic for this row.
  464.      */
  465.     protected int getBasicVariable(final int row) {
  466.         return basicRows[row];
  467.     }

  468.     /**
  469.      * Initializes the basic variable / row mapping.
  470.      * @param startColumn the column to start
  471.      */
  472.     private void initializeBasicVariables(final int startColumn) {
  473.         basicVariables = new int[getWidth() - 1];
  474.         basicRows = new int[getHeight()];

  475.         Arrays.fill(basicVariables, -1);

  476.         for (int i = startColumn; i < getWidth() - 1; i++) {
  477.             Integer row = findBasicRow(i);
  478.             if (row != null) {
  479.                 basicVariables[i] = row;
  480.                 basicRows[row] = i;
  481.             }
  482.         }
  483.     }

  484.     /**
  485.      * Returns the row in which the given column is basic.
  486.      * @param col index of the column
  487.      * @return the row that the variable is basic in, or {@code null} if the variable is not basic.
  488.      */
  489.     private Integer findBasicRow(final int col) {
  490.         Integer row = null;
  491.         for (int i = 0; i < getHeight(); i++) {
  492.             final double entry = getEntry(i, col);
  493.             if (Precision.equals(entry, 1d, maxUlps) && row == null) {
  494.                 row = i;
  495.             } else if (!Precision.equals(entry, 0d, maxUlps)) {
  496.                 return null;
  497.             }
  498.         }
  499.         return row;
  500.     }

  501.     /**
  502.      * Removes the phase 1 objective function, positive cost non-artificial variables,
  503.      * and the non-basic artificial variables from this tableau.
  504.      */
  505.     protected void dropPhase1Objective() {
  506.         if (getNumObjectiveFunctions() == 1) {
  507.             return;
  508.         }

  509.         final Set<Integer> columnsToDrop = new TreeSet<>();
  510.         columnsToDrop.add(0);

  511.         // positive cost non-artificial variables
  512.         for (int i = getNumObjectiveFunctions(); i < getArtificialVariableOffset(); i++) {
  513.             final double entry = getEntry(0, i);
  514.             if (Precision.compareTo(entry, 0d, epsilon) > 0) {
  515.                 columnsToDrop.add(i);
  516.             }
  517.         }

  518.         // non-basic artificial variables
  519.         for (int i = 0; i < getNumArtificialVariables(); i++) {
  520.             int col = i + getArtificialVariableOffset();
  521.             if (getBasicRow(col) == null) {
  522.                 columnsToDrop.add(col);
  523.             }
  524.         }

  525.         final double[][] matrix = new double[getHeight() - 1][getWidth() - columnsToDrop.size()];
  526.         for (int i = 1; i < getHeight(); i++) {
  527.             int col = 0;
  528.             for (int j = 0; j < getWidth(); j++) {
  529.                 if (!columnsToDrop.contains(j)) {
  530.                     matrix[i - 1][col++] = getEntry(i, j);
  531.                 }
  532.             }
  533.         }

  534.         // remove the columns in reverse order so the indices are correct
  535.         Integer[] drop = columnsToDrop.toArray(new Integer[0]);
  536.         for (int i = drop.length - 1; i >= 0; i--) {
  537.             columnLabels.remove((int) drop[i]);
  538.         }

  539.         this.tableau = new Array2DRowRealMatrix(matrix);
  540.         this.numArtificialVariables = 0;
  541.         // need to update the basic variable mappings as row/columns have been dropped
  542.         initializeBasicVariables(getNumObjectiveFunctions());
  543.     }

  544.     /**
  545.      * @param src the source array
  546.      * @param dest the destination array
  547.      */
  548.     private void copyArray(final double[] src, final double[] dest) {
  549.         System.arraycopy(src, 0, dest, getNumObjectiveFunctions(), src.length);
  550.     }

  551.     /**
  552.      * Returns whether the problem is at an optimal state.
  553.      * @return whether the model has been solved
  554.      */
  555.     boolean isOptimal() {
  556.         final double[] objectiveFunctionRow = getRow(0);
  557.         final int end = getRhsOffset();
  558.         for (int i = getNumObjectiveFunctions(); i < end; i++) {
  559.             final double entry = objectiveFunctionRow[i];
  560.             if (Precision.compareTo(entry, 0d, epsilon) < 0) {
  561.                 return false;
  562.             }
  563.         }
  564.         return true;
  565.     }

  566.     /**
  567.      * Get the current solution.
  568.      * @return current solution
  569.      */
  570.     protected PointValuePair getSolution() {
  571.         int negativeVarColumn = columnLabels.indexOf(NEGATIVE_VAR_COLUMN_LABEL);
  572.         Integer negativeVarBasicRow = negativeVarColumn > 0 ? getBasicRow(negativeVarColumn) : null;
  573.         double mostNegative = negativeVarBasicRow == null ? 0 : getEntry(negativeVarBasicRow, getRhsOffset());

  574.         final Set<Integer> usedBasicRows = new HashSet<>();
  575.         final double[] coefficients = new double[getOriginalNumDecisionVariables()];
  576.         for (int i = 0; i < coefficients.length; i++) {
  577.             int colIndex = columnLabels.indexOf("x" + i);
  578.             if (colIndex < 0) {
  579.                 coefficients[i] = 0;
  580.                 continue;
  581.             }
  582.             Integer basicRow = getBasicRow(colIndex);
  583.             if (basicRow != null && basicRow == 0) {
  584.                 // if the basic row is found to be the objective function row
  585.                 // set the coefficient to 0 -> this case handles unconstrained
  586.                 // variables that are still part of the objective function
  587.                 coefficients[i] = 0;
  588.             } else if (usedBasicRows.contains(basicRow)) {
  589.                 // if multiple variables can take a given value
  590.                 // then we choose the first and set the rest equal to 0
  591.                 coefficients[i] = 0 - (restrictToNonNegative ? 0 : mostNegative);
  592.             } else {
  593.                 usedBasicRows.add(basicRow);
  594.                 coefficients[i] =
  595.                     (basicRow == null ? 0 : getEntry(basicRow, getRhsOffset())) -
  596.                     (restrictToNonNegative ? 0 : mostNegative);
  597.             }
  598.             coefficients[i] = updateExponent(coefficients[i], variableExpChange[i]);
  599.         }
  600.         return new PointValuePair(coefficients, f.value(coefficients));
  601.     }

  602.     /**
  603.      * Perform the row operations of the simplex algorithm with the selected
  604.      * pivot column and row.
  605.      * @param pivotCol the pivot column
  606.      * @param pivotRow the pivot row
  607.      */
  608.     protected void performRowOperations(int pivotCol, int pivotRow) {
  609.         // set the pivot element to 1
  610.         final double pivotVal = getEntry(pivotRow, pivotCol);
  611.         divideRow(pivotRow, pivotVal);

  612.         // set the rest of the pivot column to 0
  613.         for (int i = 0; i < getHeight(); i++) {
  614.             if (i != pivotRow) {
  615.                 final double multiplier = getEntry(i, pivotCol);
  616.                 if (multiplier != 0.0) {
  617.                     subtractRow(i, pivotRow, multiplier);
  618.                 }
  619.             }
  620.         }

  621.         // update the basic variable mappings
  622.         final int previousBasicVariable = getBasicVariable(pivotRow);
  623.         basicVariables[previousBasicVariable] = -1;
  624.         basicVariables[pivotCol] = pivotRow;
  625.         basicRows[pivotRow] = pivotCol;
  626.     }

  627.     /**
  628.      * Divides one row by a given divisor.
  629.      * <p>
  630.      * After application of this operation, the following will hold:
  631.      * <pre>dividendRow = dividendRow / divisor</pre>
  632.      *
  633.      * @param dividendRowIndex index of the row
  634.      * @param divisor value of the divisor
  635.      */
  636.     protected void divideRow(final int dividendRowIndex, final double divisor) {
  637.         final double[] dividendRow = getRow(dividendRowIndex);
  638.         for (int j = 0; j < getWidth(); j++) {
  639.             dividendRow[j] /= divisor;
  640.         }
  641.     }

  642.     /**
  643.      * Subtracts a multiple of one row from another.
  644.      * <p>
  645.      * After application of this operation, the following will hold:
  646.      * <pre>minuendRow = minuendRow - multiple * subtrahendRow</pre>
  647.      *
  648.      * @param minuendRowIndex row index
  649.      * @param subtrahendRowIndex row index
  650.      * @param multiplier multiplication factor
  651.      */
  652.     protected void subtractRow(final int minuendRowIndex, final int subtrahendRowIndex, final double multiplier) {
  653.         final double[] minuendRow = getRow(minuendRowIndex);
  654.         final double[] subtrahendRow = getRow(subtrahendRowIndex);
  655.         for (int i = 0; i < getWidth(); i++) {
  656.             minuendRow[i] -= subtrahendRow[i] * multiplier;
  657.         }
  658.     }

  659.     /**
  660.      * Get the width of the tableau.
  661.      * @return width of the tableau
  662.      */
  663.     protected final int getWidth() {
  664.         return tableau.getColumnDimension();
  665.     }

  666.     /**
  667.      * Get the height of the tableau.
  668.      * @return height of the tableau
  669.      */
  670.     protected final int getHeight() {
  671.         return tableau.getRowDimension();
  672.     }

  673.     /**
  674.      * Get an entry of the tableau.
  675.      * @param row row index
  676.      * @param column column index
  677.      * @return entry at (row, column)
  678.      */
  679.     protected final double getEntry(final int row, final int column) {
  680.         return tableau.getEntry(row, column);
  681.     }

  682.     /**
  683.      * Set an entry of the tableau.
  684.      * @param row row index
  685.      * @param column column index
  686.      * @param value for the entry
  687.      */
  688.     protected final void setEntry(final int row, final int column, final double value) {
  689.         tableau.setEntry(row, column, value);
  690.     }

  691.     /**
  692.      * Get the offset of the first slack variable.
  693.      * @return offset of the first slack variable
  694.      */
  695.     protected final int getSlackVariableOffset() {
  696.         return getNumObjectiveFunctions() + numDecisionVariables;
  697.     }

  698.     /**
  699.      * Get the offset of the first artificial variable.
  700.      * @return offset of the first artificial variable
  701.      */
  702.     protected final int getArtificialVariableOffset() {
  703.         return getNumObjectiveFunctions() + numDecisionVariables + numSlackVariables;
  704.     }

  705.     /**
  706.      * Get the offset of the right hand side.
  707.      * @return offset of the right hand side
  708.      */
  709.     protected final int getRhsOffset() {
  710.         return getWidth() - 1;
  711.     }

  712.     /**
  713.      * Get the number of decision variables.
  714.      * <p>
  715.      * If variables are not restricted to positive values, this will include 1 extra decision variable to represent
  716.      * the absolute value of the most negative variable.
  717.      *
  718.      * @return number of decision variables
  719.      * @see #getOriginalNumDecisionVariables()
  720.      */
  721.     protected final int getNumDecisionVariables() {
  722.         return numDecisionVariables;
  723.     }

  724.     /**
  725.      * Get the original number of decision variables.
  726.      * @return original number of decision variables
  727.      * @see #getNumDecisionVariables()
  728.      */
  729.     protected final int getOriginalNumDecisionVariables() {
  730.         return f.getCoefficients().getDimension();
  731.     }

  732.     /**
  733.      * Get the number of slack variables.
  734.      * @return number of slack variables
  735.      */
  736.     protected final int getNumSlackVariables() {
  737.         return numSlackVariables;
  738.     }

  739.     /**
  740.      * Get the number of artificial variables.
  741.      * @return number of artificial variables
  742.      */
  743.     protected final int getNumArtificialVariables() {
  744.         return numArtificialVariables;
  745.     }

  746.     /**
  747.      * Get the row from the tableau.
  748.      * @param row the row index
  749.      * @return the reference to the underlying row data
  750.      */
  751.     protected final double[] getRow(int row) {
  752.         return tableau.getDataRef()[row];
  753.     }

  754.     /**
  755.      * Get the tableau data.
  756.      * @return tableau data
  757.      */
  758.     protected final double[][] getData() {
  759.         return tableau.getData();
  760.     }

  761.     /** {@inheritDoc} */
  762.     @Override
  763.     public boolean equals(Object other) {

  764.       if (this == other) {
  765.         return true;
  766.       }

  767.       if (other instanceof SimplexTableau) {
  768.           SimplexTableau rhs = (SimplexTableau) other;
  769.           return restrictToNonNegative  == rhs.restrictToNonNegative &&
  770.                  numDecisionVariables   == rhs.numDecisionVariables &&
  771.                  numSlackVariables      == rhs.numSlackVariables &&
  772.                  numArtificialVariables == rhs.numArtificialVariables &&
  773.                  epsilon                == rhs.epsilon &&
  774.                  maxUlps                == rhs.maxUlps &&
  775.                  f.equals(rhs.f) &&
  776.                  constraints.equals(rhs.constraints) &&
  777.                  tableau.equals(rhs.tableau);
  778.       }
  779.       return false;
  780.     }

  781.     /** {@inheritDoc} */
  782.     @Override
  783.     public int hashCode() {
  784.         return Boolean.valueOf(restrictToNonNegative).hashCode() ^
  785.                numDecisionVariables ^
  786.                numSlackVariables ^
  787.                numArtificialVariables ^
  788.                Double.valueOf(epsilon).hashCode() ^
  789.                maxUlps ^
  790.                f.hashCode() ^
  791.                constraints.hashCode() ^
  792.                tableau.hashCode();
  793.     }
  794. }