//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** @author  John Miller
 *  @version 1.0
 *  @date    Sat Aug 20 23:20:26 EDT 2011
 *  @see     LICENSE (MIT style license file).
 */

package scalation.maxima

import util.control.Breaks.{breakable, break}

import scalation.linalgebra.{MatrixD, VectorD}
import scalation.util.Error

//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** This class solves Linear Programming (LP) problems using the Revised Simplex
 *  Algorithm.  Given a constraint matrix a, constant vector b and cost vector c,
 *  find values for the solution/decision vector x that maximize the objective
 *  function f(x), while satisfying all of the constraints, i.e.,
 *
 *  maximize    f(x) = c x
 *  subject to  a x <= b, x >= 0
 *
 *  The revised algorithm has benefits over the Simplex Algorithm (less memory
 *  and reduced chance of round off errors).
 *  @see math.uc.edu/~halpern/Linear.progr.folder/Handouts.lp.02/Revisedsimplex.pdf
 *  @param a      the constraint matrix
 *  @param b      the constant/limit vector
 *  @param c      the cost/revenue vector
 *  @param x_B    the initial basis (set of indices where x_i is in the basis)
 */
class RevisedSimplex (a: MatrixD, b: VectorD, c: VectorD, var x_B: Array [Int] = null)
      extends Error
{
    private val EPSILON  = 1E-9                         // number close to zero
    private val M        = a.dim1                       // number of constraints (rows in a)
    private val N        = a.dim2                       // number of original variables (columns in a)
    private val MAX_ITER = 200 * N                      // maximum number of iterations

    if (b.dim != M) flaw ("constructor", "b.dim = " + b.dim + " != " + M)
    if (c.dim != N) flaw ("constructor", "c.dim = " + c.dim + " != " + N)

    private val b_inv    = a.selectCols (x_B).inverse   // basis of matrix-a inverted

    private var c_       = c.select (x_B) * b_inv       // adjusted cost
    private val b_       = b_inv * b                    // adjusted constants

    private var u: VectorD = null                       // ?? FIX
    private var z: VectorD = null                       // ?? FIX

    if (x_B == null) x_B = setBasis ()

    //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** There are M+N variables, N decision and M slack variables, of which,
     *  for each iteration, M are chosen for a Basic Feasible Solution (BFS).
     *  The the variables not in the basis are set to zero.  Setting j to N
     *  will start with the slack variables in the basis (only works if b >= 0).
     *  @param j  the offset to start the basis
     *  @param l  the size of the basis
     */
    def setBasis (j: Int = N, l: Int = M): Array [Int] =
    {
        val ba = Array.ofDim [Int] (l)
        for (i <- 0 until l) ba(i) = i + j
        ba
    } // setBasis

    //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Find the best variable x_l to enter the basis.
     */
    def entering (): Int = 
    {
        z = c_ * a - c
        z.argmaxPos ()
    } // entering 

    //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Find the best variable x_k to leave the basis given that x_l is entering.
     *  @param l  the variable chosen to enter the basis
     */
    def leaving (l: Int): Int =
    {
        u = b_inv * a.col(l)
        if (unbounded (u)) return -1
        var k = 0
        var r_min = Double.PositiveInfinity
        for (i <- 0 until M if u(i) > 0) {
            val r = b_(i) / u(i)
            if (r < r_min) { r_min = r; k = i}
        } // for
        k
    } // leaving

    //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Pivot by replacing x_k with x_l in the basis.  Update b_inv, b_ and c_.
     *  @param k  the leaving variable
     *  @param l  the entering variable
     */
    def pivot (k: Int, l: Int)
    {
        x_B(k) = l                // update basis (l replaces k)
        b_inv(k) /= u(k)
        b_(k)    /= u(k)
        for (i <- 0 until M if i != k) {
            b_inv(i) -= b_inv(k) * u(i)
            b_ (i)   -= b_(k) * u(i)
        } // for
        c_ -= b_inv(k) * z(l)
    } // pivot

    //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** if u <= 0., the solution is unbounded.
     *  @param u  the ?? vector FIX
     */
    def unbounded (u: VectorD): Boolean =
    {
        for (i <- 0 until u.dim if u(i) > 0.) return false
        println ("the solution is unbounded")
        true
    } // unbounded

    //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Solve a Linear Programming (LP) problem using the Revised Simplex Algorithm.
     *  Iteratively pivot until there an optimal solution is found or it is
     *  determined that the solution is unbounded.
     */
    def solve (): Tuple2 [VectorD, Double] =
    {
        var k = -1       // the leaving variable
        var l = -1       // the entering variable
        showTableau

        breakable { for (it <- 1 to MAX_ITER) {
            l = entering (); if (l == -1) break    // optimal solution found
            k = leaving (l); if (k == -1) break    // solution is unbounded
            pivot (k, l)                           // pivot: k leaves and l enters
            showTableau
        }} // for

        val x = primal
        (x, objValue (x))    // return the primal solution and the optimal value
    } // solve

    //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return the primal solution vector (x).
     */
    def primal: VectorD = b_inv * b

   //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return the dual solution vector (y).  FIX
     */
    def dual: VectorD = null

    //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return the optimal objective function value (f(x) = c x).
     *  @param x  the primal solution vector
     */
    def objValue (x: VectorD): Double = c.select(x_B) dot x

    //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Show the current revised tableau displaying the basis, b_inv, b_, c_.
     */
    def showTableau
    {
        println ("Revised Tableau")
        println ("---------------------------------------------------")
        for (i <- 0 until a.dim1) {
            println ("x" + x_B(i) + " | " + b_inv(i) + " | " + b_(i))
        } // for
        println ("c_ | " + c_)
        println ("---------------------------------------------------")
    } // showTableau

} // RevisedSimplex class


//:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** Test the Revised Simplex Algorithm class with the following maximization problem:
 *  Maximize    z = 2x_0 + 3x_1 + 4x_2
 *  Subject to      3x_0 + 2x_1 + 1x_2 + 1y_3 + 0y_4 = 10
 *                 2x_0 + 5x_1 + 3x_2 + 0y_3 + 1y_4 = 15
 *  where z is the objective variable, x are the decision variables and
 *  y are slack variables.
 */
object RevisedSimplexTest extends App
{
    // initialize matrix a, vectors b and c and the basis
    //
    val a = new MatrixD ((3, 7), 0.,  3.,  3.,  1., -1., 1., -5,
                                 1.,  3., -1.,  0., -1., 1.,  3., 
                                 1.,  2.,  0.,  1.,  2., 0., -2.)
    val c   = VectorD          (-2., -3.,  4., -3., -1., 4., -6.)
    val b   = VectorD (3., 4., 5.)
    val x_B = Array (3, 5, 0)

    val result = (new RevisedSimplex (a, b, c, x_B)).solve ()

    println ("---------------------------------------------------")
    println ("optimal x = " + result._1 + ", max cx = " + result._2)
    println ("---------------------------------------------------")

/*
    // initialize matrix a and vectors b and c
    //
    val a = new MatrixD ((2, 3), 3., 2., 1.,
                                 2., 5., 3.)
    val c   = VectorD           (2., 3., 4.)
    val b   = VectorD (10., 15.)
    val x_B = Array (1, 2)

    val result = (new RevisedSimplex (a, b, c, x_B)).solve ()

    println ("---------------------------------------------------")
    println ("optimal x = " + result._1 + ", max cx = " + result._2)
    println ("---------------------------------------------------")
*/

} // RevisedSimplexTest object