'''
The simplexLib library provides functions to load upper and lower bound, as
well as equality constraints

The library implements all steps required, i.e.,

- Translate the problem to standard format, creating auxilary variables if required
- Identify basic variables
- If basic variables are easily found, use the Big M method
- Repeat:
    - Calculate RP = relative profit for each non-basic variable
    - Calculate UL = upper limit for each basic variable
    - Replace the basic variables with the lowest UL with the non-basic with the highest RP
    - This requires a pivot-operation
    - Repeat until RP = 0 for all non-basic variable
'''

import numpy as np
# Define some flags for print...
none = 0
summary = 1
detailed = 2
nExtra = 10 # Allocate space for auxillary variables
nVar = 0
c = []
A = []
b = []
basic = [] # vector to hold basic variables
isMax = True
infty = 1e30
bigM = 100000
def extend(vec): # Routine to create extra space for auxilary variables
    for i in range(nExtra):
        vec.append(0)

def coefficients(c_in, isMaxProblem = True):
    global c,nVar
    global isMax
    isMax = isMaxProblem
    c = c_in
    nVar = len(c_in)
    if not isMaxProblem:
        for i in range(nVar):
            c[i] *= -1        
    extend(c)
    
def upperBound(b_i, A_i):
    global  nVar
    if b_i < 0:
        b_i *=-1
        for i in range(len(A_i)):
            A_i[i] *= -1
        lowerBound( b_i, A_i)  # Negative b, <= is replaced with >= when multiplying with -1
    else:
        print("Adding slack variable for upper bound constraints--> " + x(nVar) + " becomes basic variable")
        extend(A_i)        
        basic.append(nVar)
        A_i[nVar] = 1
        nVar += 1
        b.append(b_i)
        A.append(A_i)
        
def lowerBound(b_i, A_i):    
    global nVar
    if b_i < 0:
        b_i *= -1
        for i in range(len(A_i)):
            A_i[i] *= -1
        upperBound(b_i, A_i) # Negative b, >= is replaced with <= when multiplying with -1
    else:
        print("Adding surpulus variable for lower bound constraints: " + x(nVar))
        extend(A_i)
        A_i[nVar] = -1
        nVar += 1
        print("Adding basic variable with Big-M method: " + x(nVar))
        basic.append(nVar)  # Big M
        c[nVar] = -bigM
        A_i[nVar] = 1
        nVar += 1
        b.append(b_i)
        A.append(A_i)      

def equalityConstraint(b_i, A_i):
    global nVar
    if b_i < 0:
        b_i *=-1
        for i in range(len(A_i)):
            A_i[i] *= -1
    extend(A_i)
    print("Adding basic variable with Big-M method: " + x(nVar))
    basic.append(nVar)  # Big M
    c[nVar] = -bigM
    A_i[nVar] = 1
    nVar += 1
    b.append(b_i)
    A.append(A_i)

def unboundVar(vNum):
    global nVar,c,A
    print("Splitting:" + x(vNum) + " into negative and positive variable")
    c[nVar] = -c[vNum]
    for i in range(len(A)):
        A[i][nVar] = -A[i][vNum]
    nVar += 1


def SIMPLEX(printFlag=none):
    for k in range(10): # Loop until no improvement
        printProblem()
        best = -infty
        if printFlag != none:
            print("\nIteration",k+1)
        for j in range(nVar):
            test = RP(j)
            if test > best:
                best = test
                best_j = j
            if printFlag == detailed:
                print(varName(j),"RP =",test)
        if best <= 0:
            print("\nNo imporovement, final solution:")
            Z = 0
            for i in range(len(basic)):
                Z += c[basic[i]]*b[i]
                print(varName(basic[i]),"=",b[i])
            if not isMax:
                Z *=-1
            print("Z =",Z)            
            return
        best = infty
        best_i = -1
        for i in range(len(basic)):
            if A[i][best_j] == 0:
                ul = infty
            else:
                ul = b[i]/A[i][best_j]
                if ul < 0:
                    ul = infty
            if ul < best:
                best = ul
                best_i = i
            if printFlag == detailed:
                print("Basic variable",varName(basic[i]),"UL =",ul)
        if best == infty:
            print("Unbound problem")
            return     
        if printFlag != none:
            print("\nReplacing",varName(basic[best_i]), "with", varName(best_j), "as basic variable")
        pivot(best_i,best_j)
        Z = 0
        for i in range(len(basic)):
            Z += c[basic[i]]*b[i]
            print(varName(basic[i]),"=",b[i])
        if not isMax:
            Z *=-1
        print("Z (tmp)=",Z)


        
           
        
        
def RP(j):
    rp = c[j]    
    for i, basic_i in enumerate(basic):        
        rp -= c[basic_i]*A[i][j]
    return rp

def pivot(iOld,jNew):
    h = A[iOld][jNew]
    b[iOld] *= 1/h
    for j in range(nVar): # Step 1:
        A[iOld][j] *= 1/h
    for i in range(len(basic)): # Step 2:
        if i != iOld:
            h = A[i][jNew]
            b[i] -= h* b[iOld]
            for j in range(nVar):
                A[i][j] -= h * A[iOld][j]
    basic[iOld] = jNew

def varName(i):
    return "x_"+str(i+1)

def printProblem(w=6):
    print("\nProblem definition\nCoefficients=")
    printLine(c,"",w)
    print("\nConstraint matrix:")      
    for i in range(len(basic)):
          printLine(A[i],", b = " + str(round(b[i],2)),w)
          
def printLine(line,b="", w = 8):
    out=""    
    for j in range(nVar):
        if line[j] == - bigM:
            n = "-bigM"
        else:
            n = str(round(line[j],2))
        while len(n) < w:
            n = ' ' + n
        out += n
    print(out+b)

def x(n):
    return "x[" + format(n) + "] = x_" + format(n+1)