/*
 * 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.math.optimization.linear;

import java.io.IOException;
import java.io.ObjectInputStream;
import java.io.ObjectOutputStream;
import java.io.Serializable;
import java.util.ArrayList;
import java.util.Collection;
import java.util.HashSet;
import java.util.List;
import java.util.Set;

import org.apache.commons.math.linear.MatrixUtils;
import org.apache.commons.math.linear.RealMatrix;
import org.apache.commons.math.linear.Array2DRowRealMatrix;
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.MathUtils;

/**
 * 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>
 * </p>
 * @version $Revision$ $Date$
 * @since 2.0
 */
class SimplexTableau implements Serializable {

    /** Serializable version identifier. */
    private static final long serialVersionUID = -1369660067587938365L;

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

    /** Linear constraints. */
    private final Collection<LinearConstraint> constraints;

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

    /** Simple tableau. */
    protected transient RealMatrix tableau;

    /** Number of decision variables. */
    protected final int numDecisionVariables;

    /** Number of slack variables. */
    protected final int numSlackVariables;

    /** Number of artificial variables. */
    protected int numArtificialVariables;

    /** Amount of error to accept in floating point comparisons. */ 
    protected final double epsilon;
    
    /**
     * 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 in floating point comparisons
     */
    SimplexTableau(final LinearObjectiveFunction f,
                   final Collection<LinearConstraint> constraints,
                   final GoalType goalType, final boolean restrictToNonNegative,
                   final double epsilon) {
        this.f                      = f;
        this.constraints            = constraints;
        this.restrictToNonNegative  = restrictToNonNegative;
        this.epsilon                = epsilon;
        this.numDecisionVariables   = getNumVariables() + (restrictToNonNegative ? 0 : 1);
        this.numSlackVariables      = getConstraintTypeCounts(Relationship.LEQ) +
                                      getConstraintTypeCounts(Relationship.GEQ);
        this.numArtificialVariables = getConstraintTypeCounts(Relationship.EQ) +
                                      getConstraintTypeCounts(Relationship.GEQ);
        this.tableau = new Array2DRowRealMatrix(createTableau(goalType == GoalType.MAXIMIZE));
        initialize();
    }

    /**
     * Create the tableau by itself.
     * @param maximize if true, goal is to maximize the objective function
     * @return created tableau
     */
    protected double[][] createTableau(final boolean maximize) {

        // create a matrix of the correct size
        List<LinearConstraint> constraints = getNormalizedConstraints();
        int width = numDecisionVariables + numSlackVariables +
        numArtificialVariables + getNumObjectiveFunctions() + 1; // + 1 is for RHS
        int height = constraints.size() + getNumObjectiveFunctions();
        double[][] matrix = new double[height][width];

        // initialize the objective function rows
        if (getNumObjectiveFunctions() == 2) {
            matrix[0][0] = -1;
        }
        int zIndex = (getNumObjectiveFunctions() == 1) ? 0 : 1;
        matrix[zIndex][zIndex] = maximize ? 1 : -1;
        RealVector objectiveCoefficients =
            maximize ? f.getCoefficients().mapMultiply(-1) : f.getCoefficients();
            copyArray(objectiveCoefficients.getData(), matrix[zIndex], getNumObjectiveFunctions());
            matrix[zIndex][width - 1] =
                maximize ? f.getConstantTerm() : -1 * f.getConstantTerm();

                if (!restrictToNonNegative) {
                    matrix[zIndex][getSlackVariableOffset() - 1] =
                        getInvertedCoeffiecientSum(objectiveCoefficients);
                }

                // initialize the constraint rows
                int slackVar = 0;
                int artificialVar = 0;
                for (int i = 0; i < constraints.size(); i++) {
                    LinearConstraint constraint = constraints.get(i);
                    int row = getNumObjectiveFunctions() + i;

                    // decision variable coefficients
                    copyArray(constraint.getCoefficients().getData(), matrix[row], 1);

                    // x-
                    if (!restrictToNonNegative) {
                        matrix[row][getSlackVariableOffset() - 1] =
                            getInvertedCoeffiecientSum(constraint.getCoefficients());
                    }

                    // RHS
                    matrix[row][width - 1] = constraint.getValue();

                    // slack variables
                    if (constraint.getRelationship() == Relationship.LEQ) {
                        matrix[row][getSlackVariableOffset() + slackVar++] = 1;  // slack
                    } else if (constraint.getRelationship() == Relationship.GEQ) {
                        matrix[row][getSlackVariableOffset() + slackVar++] = -1; // excess
                    }

                    // artificial variables
                    if ((constraint.getRelationship() == Relationship.EQ) ||
                        (constraint.getRelationship() == Relationship.GEQ)) {
                        matrix[0][getArtificialVariableOffset() + artificialVar] = 1; 
                        matrix[row][getArtificialVariableOffset() + artificialVar++] = 1; 
                    }
                }

                return matrix;
    }

    /** Get the number of variables.
     * @return number of variables
     */
    public int getNumVariables() {
        return f.getCoefficients().getDimension();
    }

    /**
     * Get new versions of the constraints which have positive right hand sides.
     * @return new versions of the constraints
     */
    public List<LinearConstraint> getNormalizedConstraints() {
        List<LinearConstraint> normalized = new ArrayList<LinearConstraint>();
        for (LinearConstraint constraint : constraints) {
            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;
    }

    /**
     * Puts the tableau in proper form by zeroing out the artificial variables
     * in the objective function via elementary row operations.
     */
    private void initialize() {
        for (int artificialVar = 0; artificialVar < numArtificialVariables; artificialVar++) {
            int row = getBasicRow(getArtificialVariableOffset() + artificialVar);
            subtractRow(0, row, 1.0);
        }
    }

    /**
     * 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.
     */
    protected static double getInvertedCoeffiecientSum(final RealVector coefficients) {
        double sum = 0;
        for (double coefficient : coefficients.getData()) {
            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
     */
    private Integer getBasicRow(final int col) {
        Integer row = null;
        for (int i = getNumObjectiveFunctions(); i < getHeight(); i++) {
            if (MathUtils.equals(getEntry(i, col), 1.0, epsilon) && (row == null)) {
                row = i;
            } else if (!MathUtils.equals(getEntry(i, col), 0.0, epsilon)) {
                return null;
            }
        }
        return row;
    }

    /**
     * Removes the phase 1 objective function and artificial variables from this tableau.
     */
    protected void discardArtificialVariables() {
        if (numArtificialVariables == 0) {
            return;
        }
        int width = getWidth() - numArtificialVariables - 1;
        int height = getHeight() - 1;
        double[][] matrix = new double[height][width];
        for (int i = 0; i < height; i++) {
            for (int j = 0; j < width - 1; j++) {
                matrix[i][j] = getEntry(i + 1, j + 1);
            }
            matrix[i][width - 1] = getEntry(i + 1, getRhsOffset());
        }
        this.tableau = new Array2DRowRealMatrix(matrix);
        this.numArtificialVariables = 0;
    }


    /**
     * @param src the source array
     * @param dest the destination array
     * @param destPos the destination position
     */
    private void copyArray(final double[] src, final double[] dest,
                           final int destPos) {
        System.arraycopy(src, 0, dest, getNumObjectiveFunctions(), src.length); 
    }

    /**
     * Get the current solution.
     * <p>
     * {@link #solve} should be called first for this to be the optimal solution.
     * </p>
     * @return current solution
     */
    protected RealPointValuePair getSolution() {
        double[] coefficients = new double[getOriginalNumDecisionVariables()];
        Integer basicRow =
            getBasicRow(getNumObjectiveFunctions() + getOriginalNumDecisionVariables());
        double mostNegative = basicRow == null ? 0 : getEntry(basicRow, getRhsOffset());
        Set<Integer> basicRows = new HashSet<Integer>();
        for (int i = 0; i < coefficients.length; i++) {
            basicRow = getBasicRow(getNumObjectiveFunctions() + i);
            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;
            } else {
                basicRows.add(basicRow);
                coefficients[i] =
                    (basicRow == null ? 0 : getEntry(basicRow, getRhsOffset())) -
                    (restrictToNonNegative ? 0 : mostNegative);
            }
        }
        return new RealPointValuePair(coefficients, f.getValue(coefficients));
    }

    /**
     * Subtracts a multiple of one row from another.
     * <p>
     * After application of this operation, the following will hold:
     *   minuendRow = minuendRow - multiple * subtrahendRow
     * </p>
     * @param dividendRow index of the row
     * @param divisor value of the divisor
     */
    protected void divideRow(final int dividendRow, final double divisor) {
        for (int j = 0; j < getWidth(); j++) {
            tableau.setEntry(dividendRow, j, tableau.getEntry(dividendRow, j) / divisor);
        }
    }

    /**
     * Subtracts a multiple of one row from another.
     * <p>
     * After application of this operation, the following will hold:
     *   minuendRow = minuendRow - multiple * subtrahendRow
     * </p>
     * @param minuendRow row index
     * @param subtrahendRow row index
     * @param multiple multiplication factor
     */
    protected void subtractRow(final int minuendRow, final int subtrahendRow,
                               final double multiple) {
        for (int j = 0; j < getWidth(); j++) {
            tableau.setEntry(minuendRow, j, tableau.getEntry(minuendRow, j) -
                             multiple * tableau.getEntry(subtrahendRow, j));
        }
    }

    /**
     * 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.
     * </p>
     * @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 restrictToNonNegative ? numDecisionVariables : numDecisionVariables - 1;
    }

    /**
     * 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 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 == null) {
        return false;
      }

      try {

          SimplexTableau rhs = (SimplexTableau) other;
          return (restrictToNonNegative  == rhs.restrictToNonNegative) &&
                 (numDecisionVariables   == rhs.numDecisionVariables) &&
                 (numSlackVariables      == rhs.numSlackVariables) &&
                 (numArtificialVariables == rhs.numArtificialVariables) &&
                 (epsilon                == rhs.epsilon) &&
                 f.equals(rhs.f) &&
                 constraints.equals(rhs.constraints) &&
                 tableau.equals(rhs.tableau);

      } catch (ClassCastException ex) {
          // ignore exception
          return false;
      }

    }
    
    /** {@inheritDoc} */
    @Override
    public int hashCode() {
        return Boolean.valueOf(restrictToNonNegative).hashCode() ^
               numDecisionVariables ^
               numSlackVariables ^
               numArtificialVariables ^
               Double.valueOf(epsilon).hashCode() ^
               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);
    }
}
