from numpy import *
import matplotlib
from matplotlib.pyplot import *

seterr(all='print')
close('all')

sigma = 10
nprobs = 1000
eps = 1e-9
maxsize = 800

stop = 0
k = 0
ms = zeros(nprobs)
ns = zeros(nprobs)
n0s = zeros(nprobs)
iters = zeros(nprobs)
opts = zeros(nprobs)
while (k != nprobs):

    m = int(ceil(10*exp(log(maxsize/10)*random.rand())))
    n = int(ceil(10*exp(log(maxsize/10)*random.rand())))

    A = around(sigma*random.randn(m,n),0)
#    y = array(around(sigma*random.rand(1,m),0))
#    z = array(around(sigma*random.rand(1,n),0))
#    c = (matmul(y,A) - z).T
    b = array(around(sigma*abs(random.rand(m,1)),0))
    c = array(around(sigma*random.rand(n,1),0))
    
    A = -A
    AA = A
    bb = b
    cc = c
    n0 = sum(c>0)

    # 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 >= 100*(m+n):
            opt = 0
            break                      # cycling (or exponential)

    if ( (opt != 0) & (iter>0) ):
        ms[k] = m
        ns[k] = n
        n0s[k] = n0
        iters[k] = iter
        opts[k] = opt
        k = k+1
#    elif (iter == 0):                # too easy
#        #print(AA,bb,cc)
#    else:                            # cycling
#        A0 = AA[b==0,:]
#        b0 = b[b==0]
#        [m,n] = A0.shape

       
figure(1, figsize=(12, 8), dpi= 200)
plot(minimum(ms[opts>0],ns[opts>0]), iters[opts>0],'b+',label='Problems with an optimal solution')
plot(minimum(ms[opts<0],ns[opts<0]), iters[opts<0],'gx',label='Problems that are unbounded')
title('Empirical Performance of the Simplex Method')
xlabel('min(m,n)')
ylabel('number of pivots')
legend(loc='upper left', shadow=True, fontsize='x-large')

figure(2, figsize=(12, 8), dpi= 200)
loglog(ms[opts>0]+ns[opts>0], iters[opts>0],'b+',label='Problems with an optimal solution')
loglog(ms[opts<0]+ns[opts<0], iters[opts<0],'gx',label='Problems that are unbounded')
xx = arange(10,33)
yy = 2**(xx/2)
loglog(xx,yy,'r-')
ylim([1, 100000])
title('Empirical Performance of the Simplex Method')
xlabel('m+n')
ylabel('number of pivots')
legend(loc='upper left', shadow=True, fontsize='x-large')

figure(3, figsize=(12, 8), dpi= 200)
loglog(ns[opts>0], iters[opts>0],'b+',label='Problems with an optimal solution')
loglog(ns[opts<0], iters[opts<0],'gx',label='Problems that are unbounded')
title('Empirical Performance of the Simplex Method')
xlabel('n')
ylabel('number of pivots')
legend(loc='upper left', shadow=True, fontsize='x-large')

figure(4, figsize=(12, 8), dpi= 200)
loglog(1+minimum(ms[opts>0],ns[opts>0]), iters[opts>0],'b+',label='Problems with an optimal solution')
loglog(1+minimum(ms[opts<0],ns[opts<0]), iters[opts<0],'gx',label='Problems that are unbounded')
title('Empirical Performance of the Simplex Method')
xlabel('minimum of m and n')
ylabel('number of pivots')
legend(loc='upper left', shadow=True, fontsize='x-large')

logmminn = log(1+minimum(ms[opts>0],ns[opts>0]))
logiters = log(iters[opts>0])
nprobs = sum(opts>0)
X = zeros((nprobs,2))
X[:,0] = ndarray.flatten(ones((nprobs,1)))
X[:,1] = logmminn
ab = linalg.lstsq(X, logiters)[0]
a1 = exp(ab[0])
b1 = ab[1]
loglog([10, maxsize], [exp(ab[0]+ab[1]*log(10)), exp(ab[0]+ab[1]*log(maxsize))], 'b-' )

logmminn = log(1+minimum(ms[opts<0],ns[opts<0]))
logiters = log(iters[opts<0])
nprobs = sum(opts<0)
X = zeros((nprobs,2))
X[:,0] = ndarray.flatten(ones((nprobs,1)))
X[:,1] = logmminn
ab = linalg.lstsq(X, logiters)[0]
a2 = exp(ab[0])
b2 = ab[1]
loglog([10, maxsize], [exp(ab[0]+ab[1]*log(10)), exp(ab[0]+ab[1]*log(maxsize))], 'r-' )

show()