text: /*
   * Licensed to the Apache Software Foundation (ASF) under one or more
   * contributor license agreements.  See the NOTICE file distributed with
   * this work for additional information regarding copyright ownership.
   * The ASF licenses this file to You under the Apache License, Version 2.0
   * (the "License"); you may not use this file except in compliance with
   * the License.  You may obtain a copy of the License at
   *
   *      http://www.apache.org/licenses/LICENSE-2.0
   *
   * Unless required by applicable law or agreed to in writing, software
   * distributed under the License is distributed on an "AS IS" BASIS,
   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   * See the License for the specific language governing permissions and
   * limitations under the License.
   */
  package org.apache.commons.math4.optim.linear;
  
  import java.io.IOException;
  import java.io.ObjectInputStream;
  import java.io.ObjectOutputStream;
  import java.io.Serializable;
  import java.util.ArrayList;
  import java.util.Arrays;
  import java.util.Collection;
  import java.util.HashSet;
  import java.util.List;
  import java.util.Set;
  import java.util.TreeSet;
  
  import org.apache.commons.math4.exception.DimensionMismatchException;
  import org.apache.commons.math4.linear.Array2DRowRealMatrix;
  import org.apache.commons.math4.linear.MatrixUtils;
  import org.apache.commons.math4.linear.RealVector;
  import org.apache.commons.math4.optim.PointValuePair;
  import org.apache.commons.math4.optim.nonlinear.scalar.GoalType;
  import org.apache.commons.numbers.core.Precision;
  
  /**
   * A tableau for use in the Simplex method.
   *
   * <p>
   * Example:
   * <pre>
   *   W |  Z |  x1 |  x2 |  x- | s1 |  s2 |  a1 |  RHS
   * ---------------------------------------------------
   *  -1    0    0     0     0     0     0     1     0   &lt;= phase 1 objective
   *   0    1   -15   -10    0     0     0     0     0   &lt;= phase 2 objective
   *   0    0    1     0     0     1     0     0     2   &lt;= constraint 1
   *   0    0    0     1     0     0     1     0     3   &lt;= constraint 2
   *   0    0    1     1     0     0     0     1     4   &lt;= constraint 3
   * </pre>
   * W: Phase 1 objective function<br>
   * Z: Phase 2 objective function<br>
   * x1 &amp; x2: Decision variables<br>
   * x-: Extra decision variable to allow for negative values<br>
   * s1 &amp; s2: Slack/Surplus variables<br>
   * a1: Artificial variable<br>
   * RHS: Right hand side<br>
   *
   * Note on usage and safety:
   * The class is package private. It is not meant for public usage.
   * The core data structure, the tableau field, is mutated internally and
   * even reallocated when necessary.
   * Proper usage of this class is demonstrated in SimplexSolver,
   * where the class is only ever constructed in a method (never a field
   * of an object), and its lifetime, is therefore bound to a single thread (the
   * thread that's invoking the method).
   *
   * @since 2.0
   */
  class SimplexTableau implements Serializable {
  
      /** Column label for negative vars. */
      private static final String NEGATIVE_VAR_COLUMN_LABEL = "x-";
  
      /** Serializable version identifier. */
      private static final long serialVersionUID = -1369660067587938365L;
  
      /** bit mask for IEEE double exponent **/
      private static final long EXPN = 0x7ff0000000000000L;
      /** bit mask for IEEE double mantissa and sign **/
      private static final long FRAC = 0x800fffffffffffffL;
      /** max IEEE exponent is 2047 **/
      private static final int MAX_IEEE_EXP = 2047;
      /** min IEEE exponent is 0 **/
      private static final int MIN_IEEE_EXP = 0;
      /** IEEE exponent is kept in an offset form, 1023 is zero **/
      private static final int OFFSET_IEEE_EXP = 1023;
      /** double exponent shift per IEEE standard **/
      private static final int IEEE_EXPONENT_SHIFT = 52;
  
      /** Linear objective function. */
      private final LinearObjectiveFunction f;
  
      /** Linear constraints. */
      private final List<LinearConstraint> constraints;
  
      /** Whether to restrict the variables to non-negative values. */
      private final boolean restrictToNonNegative;
  
      /** The variables each column represents */
      private final List<String> columnLabels = new ArrayList<>();
  
      /** Simple tableau. */
      private transient Array2DRowRealMatrix tableau;
  
      /** Number of decision variables. */
      private final int numDecisionVariables;
  
      /** Number of slack variables. */
      private final int numSlackVariables;
  
      /** Number of artificial variables. */
      private int numArtificialVariables;
  
      /** 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;
  
      /** Maps basic variables to row they are basic in. */
      private int[] basicVariables;
  
      /** Maps rows to their corresponding basic variables. */
      private int[] basicRows;
  
      /** changes in floating point exponent to scale the input */
      private int[] variableExpChange;
  
      /**
       * Builds a tableau for a linear problem.
       *
       * @param f Linear objective function.
       * @param constraints Linear constraints.
       * @param goalType 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.
       * @throws DimensionMismatchException if the dimension of the constraints does not match the
       *   dimension of the objective function
       */
      SimplexTableau(final LinearObjectiveFunction f,
                     final Collection<LinearConstraint> constraints,
                     final GoalType goalType,
                     final boolean restrictToNonNegative,
                     final double epsilon) {
          this(f, constraints, goalType, restrictToNonNegative, epsilon, SimplexSolver.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
       * @throws DimensionMismatchException if the dimension of the constraints does not match the
       *   dimension of the objective function
       */
      SimplexTableau(final LinearObjectiveFunction f,
                     final Collection<LinearConstraint> constraints,
                     final GoalType goalType,
                     final boolean restrictToNonNegative,
                     final double epsilon,
                     final int maxUlps) throws DimensionMismatchException {
          checkDimensions(f, constraints);
          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) +
                                        getConstraintTypeCounts(Relationship.GEQ);
          this.numArtificialVariables = getConstraintTypeCounts(Relationship.EQ) +
                                        getConstraintTypeCounts(Relationship.GEQ);
          this.tableau = createTableau(goalType == GoalType.MAXIMIZE);
          // initialize the basic variables for phase 1:
          //   we know that only slack or artificial variables can be basic
          initializeBasicVariables(getSlackVariableOffset());
          initializeColumnLabels();
      }
  
      /**
       * Checks that the dimensions of the objective function and the constraints match.
       * @param objectiveFunction the objective function
       * @param c the set of constraints
       * @throws DimensionMismatchException if the constraint dimensions do not match with the
       *   dimension of the objective function
       */
      private void checkDimensions(final LinearObjectiveFunction objectiveFunction,
                                   final Collection<LinearConstraint> c) {
          final int dimension = objectiveFunction.getCoefficients().getDimension();
          for (final LinearConstraint constraint : c) {
              final int constraintDimension = constraint.getCoefficients().getDimension();
              if (constraintDimension != dimension) {
                  throw new DimensionMismatchException(constraintDimension, dimension);
              }
          }
      }
      /**
       * Initialize the labels for the columns.
       */
      protected void initializeColumnLabels() {
        if (getNumObjectiveFunctions() == 2) {
          columnLabels.add("W");
        }
        columnLabels.add("Z");
        for (int i = 0; i < getOriginalNumDecisionVariables(); i++) {
          columnLabels.add("x" + i);
        }
        if (!restrictToNonNegative) {
          columnLabels.add(NEGATIVE_VAR_COLUMN_LABEL);
        }
        for (int i = 0; i < getNumSlackVariables(); i++) {
          columnLabels.add("s" + i);
        }
        for (int i = 0; i < getNumArtificialVariables(); i++) {
          columnLabels.add("a" + i);
        }
        columnLabels.add("RHS");
      }
  
      /**
       * Create the tableau by itself.
       * @param maximize if true, goal is to maximize the objective function
       * @return created tableau
       */
      protected Array2DRowRealMatrix createTableau(final boolean maximize) {
  
          // create a matrix of the correct size
          int width = numDecisionVariables + numSlackVariables +
          numArtificialVariables + getNumObjectiveFunctions() + 1; // + 1 is for RHS
          int height = constraints.size() + getNumObjectiveFunctions();
          Array2DRowRealMatrix matrix = new Array2DRowRealMatrix(height, width);
  
          // initialize the objective function rows
          if (getNumObjectiveFunctions() == 2) {
              matrix.setEntry(0, 0, -1);
          }
  
          int zIndex = (getNumObjectiveFunctions() == 1) ? 0 : 1;
          matrix.setEntry(zIndex, zIndex, maximize ? 1 : -1);
  
          double[][] scaled = new double[constraints.size() + 1][];
  
          RealVector objectiveCoefficients = maximize ? f.getCoefficients().mapMultiply(-1) : f.getCoefficients();
          scaled[0] = objectiveCoefficients.toArray();
          double[] scaledRhs = new double[constraints.size() + 1];
          double value = maximize ? f.getConstantTerm() : -1 * f.getConstantTerm();
          scaledRhs[0] = value;
  
          for (int i = 0; i < constraints.size(); i++) {
              LinearConstraint constraint = constraints.get(i);
              scaled[i + 1] = constraint.getCoefficients().toArray();
              scaledRhs[i + 1] = constraint.getValue();
          }
          variableExpChange = new int[scaled[0].length];
  
          scale(scaled, scaledRhs);
  
          copyArray(scaled[0], matrix.getDataRef()[zIndex]);
          matrix.setEntry(zIndex, width - 1, scaledRhs[0]);
  
          if (!restrictToNonNegative) {
              matrix.setEntry(zIndex, getSlackVariableOffset() - 1,
                              getInvertedCoefficientSum(scaled[0]));
          }
  
          // initialize the constraint rows
          int slackVar = 0;
          int artificialVar = 0;
          for (int i = 0; i < constraints.size(); i++) {
              final LinearConstraint constraint = constraints.get(i);
              final int row = getNumObjectiveFunctions() + i;
  
              // decision variable coefficients
              copyArray(scaled[i + 1], matrix.getDataRef()[row]);
  
              // x-
              if (!restrictToNonNegative) {
                  matrix.setEntry(row, getSlackVariableOffset() - 1,
                                  getInvertedCoefficientSum(scaled[i + 1]));
              }
  
              // RHS
              matrix.setEntry(row, width - 1, scaledRhs[i + 1]);
  
              // slack variables
              if (constraint.getRelationship() == Relationship.LEQ) {
                  matrix.setEntry(row, getSlackVariableOffset() + slackVar++, 1);  // slack
              } else if (constraint.getRelationship() == Relationship.GEQ) {
                  matrix.setEntry(row, getSlackVariableOffset() + slackVar++, -1); // excess
              }
  
              // artificial variables
              if ((constraint.getRelationship() == Relationship.EQ) ||
                  (constraint.getRelationship() == Relationship.GEQ)) {
                  matrix.setEntry(0, getArtificialVariableOffset() + artificialVar, 1);
                  matrix.setEntry(row, getArtificialVariableOffset() + artificialVar++, 1);
                  matrix.setRowVector(0, matrix.getRowVector(0).subtract(matrix.getRowVector(row)));
              }
          }
  
          return matrix;
      }
  
      /** We scale the constants in the equations and objective, which means we try
       * to get the IEEE double exponent as close to zero (1023) as possible, which makes the
       * constants closer to 1.
       * We use exponent shifts instead of division because that introduces no bit errors.
       *
       * @param scaled coefficients before scaling
       * @param scaledRhs right hand side before scaling
       */
      private void scale(double[][] scaled, double[] scaledRhs) {
          /*
              first transform across:
              c0 x0 + c1 x1 + ... + cn xn = vn ==> (2^expChange) * (c0 x0 + c1 x1 + ... + cn xn = vn)
  
              expChange will be negative if the constants are larger than 1,
              it'll be positive if the constants are less than 1.
          */
          for (int i = 0; i < scaled.length; i++) {
              int minExp = MAX_IEEE_EXP + 1;
              int maxExp = MIN_IEEE_EXP - 1;
              for (double d: scaled[i]) {
                  if (d != 0) {
                      int e = exponent(d);
                      if (e < minExp) {
                          minExp = e;
                      }
                      if (e > maxExp) {
                          maxExp = e;
                      }
                  }
              }
              if (scaledRhs[i] != 0) {
                  final int e = exponent(scaledRhs[i]);
                  if (e < minExp) {
                      minExp = e;
                  }
                  if (e > maxExp) {
                      maxExp = e;
                  }
              }
              final int expChange = computeExpChange(minExp, maxExp);
              if (expChange != 0) {
                  scaledRhs[i] = updateExponent(scaledRhs[i], expChange);
                  updateExponent(scaled[i], expChange);
              }
          }
  
          /*
              second, transform down the columns. this is like defining a new variable for that column
              that is yi = xi * (2^expChange)
              After solving for yi, we compute xi by shifting again. See getSolution()
           */
          for (int i = 0; i < variableExpChange.length; i++) {
              int minExp = MAX_IEEE_EXP + 1;
              int maxExp = MIN_IEEE_EXP - 1;
  
              for (double[] coefficients : scaled) {
                  final double d = coefficients[i];
                  if (d != 0) {
                      int e = exponent(d);
                      if (e < minExp) {
                          minExp = e;
                      }
                      if (e > maxExp) {
                          maxExp = e;
                      }
                  }
              }
              final int expChange = computeExpChange(minExp, maxExp);
              variableExpChange[i] = expChange;
              if (expChange != 0) {
                  for (double[] coefficients : scaled) {
                       coefficients[i] = updateExponent(coefficients[i], expChange);
                  }
              }
          }
      }
  
      /**
       * Given the minimum and maximum value of the exponent of two {@code double}
       * values, pick a change in exponent to bring those values closer to 1.
       *
       * @param minExp Smallest exponent.
       * @param maxExp Largest exponent.
       * @return the new exponent.
       */
      private int computeExpChange(int minExp, int maxExp) {
          int expChange = 0;
          if (minExp <= MAX_IEEE_EXP &&
              minExp > OFFSET_IEEE_EXP) {
              expChange = OFFSET_IEEE_EXP - minExp;
          } else if (maxExp >= MIN_IEEE_EXP &&
                     maxExp < OFFSET_IEEE_EXP) {
              expChange = OFFSET_IEEE_EXP - maxExp;
          }
          return expChange;
      }
  
      /**
       * Changes the exponent of every member of the array by the given amount.
       *
       * @param dar array of doubles to change
       * @param exp exponent value to change
       */
      private static void updateExponent(double[] dar, int exp) {
          for (int i = 0; i < dar.length; i++) {
              dar[i] = updateExponent(dar[i], exp);
          }
      }
  
      /**
       * Extract the exponent of a {@code double}.
       *
       * @param d value to extract the exponent from
       * @return the IEEE exponent in the EXPN bits, as an integer
       */
      private static int exponent(double d) {
          final long bits = Double.doubleToLongBits(d);
          return (int) ((bits & EXPN) >>> IEEE_EXPONENT_SHIFT);
      }
  
      /**
       * Changes the exponent of a number by the given amount.
       *
       * @param d value to change
       * @param exp exponent to add to the existing exponent (may be negative)
       * @return a double with the same sign/mantissa bits as d, but exponent changed by exp
       */
      private static double updateExponent(double d, int exp) {
          if (d == 0 ||
              exp == 0) {
              return d;
          }
          final long bits = Double.doubleToLongBits(d);
          return Double.longBitsToDouble((bits & FRAC) | ((((bits & EXPN) >>> IEEE_EXPONENT_SHIFT) + exp) << IEEE_EXPONENT_SHIFT));
      }
  
      /**
       * Get new versions of the constraints which have positive right hand sides.
       * @param originalConstraints original (not normalized) constraints
       * @return new versions of the constraints
       */
      public List<LinearConstraint> normalizeConstraints(Collection<LinearConstraint> originalConstraints) {
          final List<LinearConstraint> normalized = new ArrayList<>(originalConstraints.size());
          for (LinearConstraint constraint : originalConstraints) {
              normalized.add(normalize(constraint));
          }
          return normalized;
      }
  
      /**
       * Get a new equation equivalent to this one with a positive right hand side.
       * @param constraint reference constraint
       * @return new equation
       */
      private LinearConstraint normalize(final LinearConstraint constraint) {
          if (constraint.getValue() < 0) {
              return new LinearConstraint(constraint.getCoefficients().mapMultiply(-1),
                                          constraint.getRelationship().oppositeRelationship(),
                                          -1 * constraint.getValue());
          }
          return new LinearConstraint(constraint.getCoefficients(),
                                      constraint.getRelationship(), constraint.getValue());
      }
  
      /**
       * Get the number of objective functions in this tableau.
       * @return 2 for Phase 1.  1 for Phase 2.
       */
      protected final int getNumObjectiveFunctions() {
          return this.numArtificialVariables > 0 ? 2 : 1;
      }
  
      /**
       * Get a count of constraints corresponding to a specified relationship.
       * @param relationship relationship to count
       * @return number of constraint with the specified relationship
       */
      private int getConstraintTypeCounts(final Relationship relationship) {
          int count = 0;
          for (final LinearConstraint constraint : constraints) {
              if (constraint.getRelationship() == relationship) {
                  ++count;
              }
          }
          return count;
      }
  
      /**
       * Get the -1 times the sum of all coefficients in the given array.
       * @param coefficients coefficients to sum
       * @return the -1 times the sum of all coefficients in the given array.
       */
      private static double getInvertedCoefficientSum(final double[] coefficients) {
          double sum = 0;
          for (double coefficient : coefficients) {
              sum -= coefficient;
          }
          return sum;
      }
  
      /**
       * Checks whether the given column is basic.
       * @param col index of the column to check
       * @return the row that the variable is basic in.  null if the column is not basic
       */
      protected Integer getBasicRow(final int col) {
          final int row = basicVariables[col];
          return row == -1 ? null : row;
      }
  
      /**
       * Returns the variable that is basic in this row.
       * @param row the index of the row to check
       * @return the variable that is basic for this row.
       */
      protected int getBasicVariable(final int row) {
          return basicRows[row];
      }
  
      /**
       * Initializes the basic variable / row mapping.
       * @param startColumn the column to start
       */
      private void initializeBasicVariables(final int startColumn) {
          basicVariables = new int[getWidth() - 1];
          basicRows = new int[getHeight()];
  
          Arrays.fill(basicVariables, -1);
  
          for (int i = startColumn; i < getWidth() - 1; i++) {
              Integer row = findBasicRow(i);
              if (row != null) {
                  basicVariables[i] = row;
                  basicRows[row] = i;
              }
          }
      }
  
      /**
       * Returns the row in which the given column is basic.
       * @param col index of the column
       * @return the row that the variable is basic in, or {@code null} if the variable is not basic.
       */
      private Integer findBasicRow(final int col) {
          Integer row = null;
          for (int i = 0; i < getHeight(); i++) {
              final double entry = getEntry(i, col);
              if (Precision.equals(entry, 1d, maxUlps) && (row == null)) {
                  row = i;
              } else if (!Precision.equals(entry, 0d, maxUlps)) {
                  return null;
              }
          }
          return row;
      }
  
      /**
       * Removes the phase 1 objective function, positive cost non-artificial variables,
       * and the non-basic artificial variables from this tableau.
       */
      protected void dropPhase1Objective() {
          if (getNumObjectiveFunctions() == 1) {
              return;
          }
  
          final Set<Integer> columnsToDrop = new TreeSet<>();
          columnsToDrop.add(0);
  
          // positive cost non-artificial variables
          for (int i = getNumObjectiveFunctions(); i < getArtificialVariableOffset(); i++) {
              final double entry = getEntry(0, i);
              if (Precision.compareTo(entry, 0d, epsilon) > 0) {
                  columnsToDrop.add(i);
              }
          }
  
          // non-basic artificial variables
          for (int i = 0; i < getNumArtificialVariables(); i++) {
              int col = i + getArtificialVariableOffset();
              if (getBasicRow(col) == null) {
                  columnsToDrop.add(col);
              }
          }
  
          final double[][] matrix = new double[getHeight() - 1][getWidth() - columnsToDrop.size()];
          for (int i = 1; i < getHeight(); i++) {
              int col = 0;
              for (int j = 0; j < getWidth(); j++) {
                  if (!columnsToDrop.contains(j)) {
                      matrix[i - 1][col++] = getEntry(i, j);
                  }
              }
          }
  
          // remove the columns in reverse order so the indices are correct
          Integer[] drop = columnsToDrop.toArray(new Integer[columnsToDrop.size()]);
          for (int i = drop.length - 1; i >= 0; i--) {
              columnLabels.remove((int) drop[i]);
          }
  
          this.tableau = new Array2DRowRealMatrix(matrix);
          this.numArtificialVariables = 0;
          // need to update the basic variable mappings as row/columns have been dropped
          initializeBasicVariables(getNumObjectiveFunctions());
      }
  
      /**
       * @param src the source array
       * @param dest the destination array
       */
      private void copyArray(final double[] src, final double[] dest) {
          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() {
          final double[] objectiveFunctionRow = getRow(0);
          final int end = getRhsOffset();
          for (int i = getNumObjectiveFunctions(); i < end; i++) {
              final double entry = objectiveFunctionRow[i];
              if (Precision.compareTo(entry, 0d, epsilon) < 0) {
                  return false;
              }
          }
          return true;
      }
  
      /**
       * Get the current solution.
       * @return current solution
       */
      protected PointValuePair getSolution() {
          int negativeVarColumn = columnLabels.indexOf(NEGATIVE_VAR_COLUMN_LABEL);
          Integer negativeVarBasicRow = negativeVarColumn > 0 ? getBasicRow(negativeVarColumn) : null;
          double mostNegative = negativeVarBasicRow == null ? 0 : getEntry(negativeVarBasicRow, getRhsOffset());
  
          final Set<Integer> usedBasicRows = new HashSet<>();
          final double[] coefficients = new double[getOriginalNumDecisionVariables()];
          for (int i = 0; i < coefficients.length; i++) {
              int colIndex = columnLabels.indexOf("x" + i);
              if (colIndex < 0) {
                  coefficients[i] = 0;
                  continue;
              }
              Integer basicRow = getBasicRow(colIndex);
              if (basicRow != null && basicRow == 0) {
                  // if the basic row is found to be the objective function row
                  // set the coefficient to 0 -> this case handles unconstrained
                  // variables that are still part of the objective function
                  coefficients[i] = 0;
              } else if (usedBasicRows.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 - (restrictToNonNegative ? 0 : mostNegative);
              } else {
                  usedBasicRows.add(basicRow);
                  coefficients[i] =
                      (basicRow == null ? 0 : getEntry(basicRow, getRhsOffset())) -
                      (restrictToNonNegative ? 0 : mostNegative);
              }
              coefficients[i] = updateExponent(coefficients[i], variableExpChange[i]);
          }
          return new PointValuePair(coefficients, f.value(coefficients));
      }
  
      /**
       * Perform the row operations of the simplex algorithm with the selected
       * pivot column and row.
       * @param pivotCol the pivot column
       * @param pivotRow the pivot row
       */
      protected void performRowOperations(int pivotCol, int pivotRow) {
          // set the pivot element to 1
          final double pivotVal = getEntry(pivotRow, pivotCol);
          divideRow(pivotRow, pivotVal);
  
          // set the rest of the pivot column to 0
          for (int i = 0; i < getHeight(); i++) {
              if (i != pivotRow) {
                  final double multiplier = getEntry(i, pivotCol);
                  if (multiplier != 0.0) {
                      subtractRow(i, pivotRow, multiplier);
                  }
              }
          }
  
          // update the basic variable mappings
          final int previousBasicVariable = getBasicVariable(pivotRow);
          basicVariables[previousBasicVariable] = -1;
          basicVariables[pivotCol] = pivotRow;
          basicRows[pivotRow] = pivotCol;
      }
  
      /**
       * Divides one row by a given divisor.
       * <p>
       * After application of this operation, the following will hold:
       * <pre>dividendRow = dividendRow / divisor</pre>
       *
       * @param dividendRowIndex index of the row
       * @param divisor value of the divisor
       */
      protected void divideRow(final int dividendRowIndex, final double divisor) {
          final double[] dividendRow = getRow(dividendRowIndex);
          for (int j = 0; j < getWidth(); j++) {
              dividendRow[j] /= divisor;
          }
      }
  
      /**
       * Subtracts a multiple of one row from another.
       * <p>
       * After application of this operation, the following will hold:
       * <pre>minuendRow = minuendRow - multiple * subtrahendRow</pre>
       *
       * @param minuendRowIndex row index
       * @param subtrahendRowIndex row index
       * @param multiplier multiplication factor
       */
      protected void subtractRow(final int minuendRowIndex, final int subtrahendRowIndex, final double multiplier) {
          final double[] minuendRow = getRow(minuendRowIndex);
          final double[] subtrahendRow = getRow(subtrahendRowIndex);
          for (int i = 0; i < getWidth(); i++) {
              minuendRow[i] -= subtrahendRow[i] * multiplier;
          }
      }
  
      /**
       * Get the width of the tableau.
       * @return width of the tableau
       */
      protected final int getWidth() {
          return tableau.getColumnDimension();
      }
  
      /**
       * Get the height of the tableau.
       * @return height of the tableau
       */
      protected final int getHeight() {
          return tableau.getRowDimension();
      }
  
      /**
       * Get an entry of the tableau.
       * @param row row index
       * @param column column index
       * @return entry at (row, column)
       */
      protected final double getEntry(final int row, final int column) {
          return tableau.getEntry(row, column);
      }
  
      /**
       * Set an entry of the tableau.
       * @param row row index
       * @param column column index
       * @param value for the entry
       */
      protected final void setEntry(final int row, final int column, final double value) {
          tableau.setEntry(row, column, value);
      }
  
      /**
       * Get the offset of the first slack variable.
       * @return offset of the first slack variable
       */
      protected final int getSlackVariableOffset() {
          return getNumObjectiveFunctions() + numDecisionVariables;
      }
  
      /**
       * Get the offset of the first artificial variable.
       * @return offset of the first artificial variable
       */
      protected final int getArtificialVariableOffset() {
          return getNumObjectiveFunctions() + numDecisionVariables + numSlackVariables;
      }
  
      /**
       * Get the offset of the right hand side.
       * @return offset of the right hand side
       */
      protected final int getRhsOffset() {
          return getWidth() - 1;
      }
  
      /**
       * Get the number of decision variables.
       * <p>
       * If variables are not restricted to positive values, this will include 1 extra decision variable to represent
       * the absolute value of the most negative variable.
       *
       * @return number of decision variables
       * @see #getOriginalNumDecisionVariables()
       */
      protected final int getNumDecisionVariables() {
          return numDecisionVariables;
      }
  
      /**
       * Get the original number of decision variables.
       * @return original number of decision variables
       * @see #getNumDecisionVariables()
       */
      protected final int getOriginalNumDecisionVariables() {
          return f.getCoefficients().getDimension();
      }
  
      /**
       * Get the number of slack variables.
       * @return number of slack variables
       */
      protected final int getNumSlackVariables() {
          return numSlackVariables;
      }
  
      /**
       * Get the number of artificial variables.
       * @return number of artificial variables
       */
      protected final int getNumArtificialVariables() {
          return numArtificialVariables;
      }
  
      /**
       * Get the row from the tableau.
       * @param row the row index
       * @return the reference to the underlying row data
       */
      protected final double[] getRow(int row) {
          return tableau.getDataRef()[row];
      }
  
      /**
       * Get the tableau data.
       * @return tableau data
       */
      protected final double[][] getData() {
          return tableau.getData();
      }
  
      /** {@inheritDoc} */
      @Override
      public boolean equals(Object other) {
  
        if (this == other) {
          return true;
        }
  
        if (other instanceof SimplexTableau) {
            SimplexTableau rhs = (SimplexTableau) other;
            return (restrictToNonNegative  == rhs.restrictToNonNegative) &&
                   (numDecisionVariables   == rhs.numDecisionVariables) &&
                   (numSlackVariables      == rhs.numSlackVariables) &&
                   (numArtificialVariables == rhs.numArtificialVariables) &&
                   (epsilon                == rhs.epsilon) &&
                   (maxUlps                == rhs.maxUlps) &&
                   f.equals(rhs.f) &&
                   constraints.equals(rhs.constraints) &&
                   tableau.equals(rhs.tableau);
        }
        return false;
      }
  
      /** {@inheritDoc} */
      @Override
      public int hashCode() {
          return Boolean.valueOf(restrictToNonNegative).hashCode() ^
                 numDecisionVariables ^
                 numSlackVariables ^
                 numArtificialVariables ^
                 Double.valueOf(epsilon).hashCode() ^
                 maxUlps ^
                 f.hashCode() ^
                 constraints.hashCode() ^
                 tableau.hashCode();
      }
  
      /**
       * Serialize the instance.
       * @param oos stream where object should be written
       * @throws IOException if object cannot be written to stream
       */
      private void writeObject(ObjectOutputStream oos)
          throws IOException {
          oos.defaultWriteObject();
          MatrixUtils.serializeRealMatrix(tableau, oos);
      }
  
      /**
       * Deserialize the instance.
       * @param ois stream from which the object should be read
       * @throws ClassNotFoundException if a class in the stream cannot be found
       * @throws IOException if object cannot be read from the stream
       */
      private void readObject(ObjectInputStream ois)
        throws ClassNotFoundException, IOException {
          ois.defaultReadObject();
          MatrixUtils.deserializeRealMatrix(this, "tableau", ois);
      }
  }
  