# this is an adaptation of the primalSimplex.py provided in www ref of the course textbook by prof. R. J. Vanderbei

from numpy import *

nprobs = 1000000
eps = 1e-9

k = 1
iters = zeros(nprobs)
while (k != nprobs):
    
    if k%100000 == 0 :
        print(k, 'instance generated')
    
    m = 2
    n = 4

    A = around(100*random.randn(m,n),2)
    b = zeros((m,1))
    c = around(100*random.rand(n,1),2)
    
    A = -A

    # here's the primal simplex algorithm
    #
    # Problem:    max c^T x
    #             s.t. A  x <= b
    #                     x >= 0
    #
    nonbasics = arange(n)
    basics = n + arange(m)
    
    iter = 0
    opt = 1
    while ( (max(c) > eps) ):
        col = argmax(c)                 # find entering variable
        Acol = A[:,col].reshape(m,1)    # the associate entering column
        tmp = -Acol/(b+eps)             # vector of ratios
        row = argmax(tmp)               # the leaving variable

        if ( sum( Acol<-eps ) == 0 ): 
            opt = -1                    # unbounded
            break

        j = nonbasics[col]              # the entering var's subscript
        i = basics[row]                 # the leaving var's subscript

        # update A matrix.  See section 5.4.
        Arow = A[row,:]                 # the row in A of the entering variable
        a = A[row,col]                  # the pivot element
        A = A - Acol*Arow/a             # the out-of-row/out-of-col update formula
        A[row,:] = -Arow/a              # update formula for the row
        A[:,col] = Acol.reshape(1,m)/a  # update formula for the col
        A[row,col] = 1/a                # update formula for the pivot element

        # update the right-hand side
        brow = b[row,0]
        b = b - brow*Acol/a
        b[row] = -brow/a

        # update the objective function
        ccol = c[col,0]
        c = c - ccol*(Arow.reshape(n,1))/a
        c[col] = ccol/a

        # swapping variables x_j and x_i position in the dictionary
        basics[row] = j         
        nonbasics[col] = i
        iter = iter+1
        if iter >= 15:
            opt = 0
            print("found a cycling example");
            break                      # cycling (or exponential)
    
    if ( (opt != 0) & (iter>0) ):
        iters[k] = iter
        k = k+1


print( max(iters) )