# All indices (of vectors, matrices) are assumed to start at 0

def from_list(b,L):
    c = vector(QQ, len(L))
    j = 0
    for i in L:
        c[j]=b[i]
        j = j+1
    return c

def ratios(z,y):
    r = vector(QQ,len(z))
    for i in range(len(z)):
        if y[i]<0 :
            r[i]= (-z[i] / y[i])
        else:
            r[i] = -1 

    return r

# Compute a maximal independent subset of the rows of A that include B.
# Note that A is assumed to have full column rank.

def ind_rows(A,B):
	L = deepcopy(B)
	for i in range(0,A.nrows()):
		L.append(i)
		if A.matrix_from_rows(L).rank() > A.matrix_from_rows(B).rank():
			B.append(i)
		else:
			L.remove(i)
	return B

def step_3(A,c,B,i):              

    z = A.matrix_from_rows(B).solve_left(c)
    y = A.matrix_from_rows(B).solve_left(-A.row(i))

    r = ratios(z,y)

    l = infinity
    j = -1
    for k in range (A.ncols()):

       if (r[k] >= 0):
           if (r[k] < l):
               l = r[k]
               j = k
    
    if (j < 0):
        return [B,False]

    B[j] = i

    return [B,True]

# input: matrix A, initial roof B, right hand side b, objective c
# output: either feasible solution and a certificate for its optimality w.r.t. to c
#         or an empty vector and an infeasibility certificate

def phase_II(A,B,b,c):
	while True:
		vertex = A.matrix_from_rows(B).solve_right(from_list(b,B))
		# sort the current roof (Least entering and exiting rule)
		B.sort()
		diff = A*vertex - b
		i_enter = -1
		for i in range(0, len(diff)):
			if diff[i] > 0:
				i_enter = i
				break
		if i_enter < 0:
			# vertex is feasible
			l = A.matrix_from_rows(B).solve_left(c)
			l_ext = vector(QQ,A.nrows())
			for j in range(0,len(B)):
				l_ext[B[j]] = l[j]
			return [vertex,l_ext,B]
		else:
			[B,feasible] = step_3(A,c,B,i_enter)
			if not(feasible):
				B.append(i_enter)
				#find infeasibility certificate
				Ker = A.matrix_from_rows(B).left_kernel()
				if len(Ker.basis()) != 1:
					print "Error in computing infeasibility certificate"
				l = Ker.basis_matrix()[0]

				norm = abs(l*from_list(b,B))
				l_ext = vector(QQ,A.nrows())
				for j in range(0,len(B)):
					l_ext[B[j]] = l[j]/norm
				B.remove(i_enter)
				return [vector(QQ,0),l_ext,B]

# input: matrix A, objective c
# output: either an initial roof
#         or indication that the LP is infeasible or unbounded.

def phase_I(A,c):
	if c.is_zero():
		return [ind_rows(A,[]),True,vector(QQ,0)]

	n = A.ncols()
	m = A.nrows()
	A_ext = A.stack (matrix(c))
	b_ext = vector(QQ,m+1)
	b_ext[m] = 1

	B = ind_rows(A_ext,[m])
	[x,l,B] = phase_II(A_ext,B,b_ext,c)

	# is c^Tx <= 1 in the optimal roof? Or is x = 0 optimal solution?
	if c*x > 0:
		return [[],False,x]
	else:
		return [B,True,vector(QQ,0)]

def gauss2(A):               
    U = matrix(QQ,A.nrows(),A.nrows())
    for i in range(A.nrows()):
        U[i,i]=1      
        
    U1 = deepcopy(U)

    k = 1
    for j in range(1, A.ncols()+1):
        i=k        
        while(i <= A.nrows() and A[i-1,j-1] == 0):
            i = i+1
        if (i <= A.nrows() ):
            U.swap_rows(k-1,i-1)
            U1.swap_columns(k-1,i-1)
            A.swap_rows(k-1,i-1)

            U.rescale_row(k-1, 1/A[k-1,j-1])
            U1.rescale_col(k-1, A[k-1,j-1])
            A.rescale_row(k-1, 1/A[k-1,j-1])           

                
            for i in range(k+1,A.nrows()+1):
                U1.add_multiple_of_column(k-1,i-1,A[i-1,j-1])
                U.add_multiple_of_row(i-1,k-1,-A[i-1,j-1])
                A.add_multiple_of_row(i-1,k-1,-A[i-1,j-1])


            for i in range(1,k):
                U1.add_multiple_of_column(k-1,i-1,A[i-1,j-1])
                U.add_multiple_of_row(i-1,k-1,-A[i-1,j-1])
                A.add_multiple_of_row(i-1,k-1,-A[i-1,j-1])
                
            k = k+1

    return (U,U1,k-1)

# input: matrix A, objective c
# output: a matrix of full column rank and the corresponding transformed objective
#         indication if the LP is unbounded (Attention: No sufficient condition for the boundedness of the LP)

def preprocessing(A,c):

	Aprime = deepcopy(A.transpose())
	(U,U1,k) = gauss2(Aprime)
	A = A*U.transpose()
	c = c*U.transpose()	

	B = [i for i in range(0,k)]
	A = A.matrix_from_columns(B)

	bounded = True
	for i in range(k, len(c)):
		if c[i] != 0:
			bounded = False
	c = from_list(c,B)

	return [A,c,U,bounded]

def simplex(M,v,w):

	A = deepcopy(M)
	b = deepcopy(v)
	c = deepcopy(w)

	[A,c,U,bounded] = preprocessing(A,c)

	[B,found,d] = phase_I(A,c)
	if found:
		[x,l,B] = phase_II(A,B,b,c)

		if len(x) == 0:
			print "LP is infeasible"
			print "Check:"
			print l*M == 0 and l*v == -1 
			return [x,l]
		else:
			x_ext = vector(QQ,M.ncols())
			for i in range(0,len(x)):
				x_ext[i] = x[i]
			if bounded:
				print "LP is feasible and bounded"
				print "Check:"
				print M*U.transpose()*x_ext <= v and l*M == w
				return [U.transpose()*x_ext,l]
			else:
				# compute unbounded certificate
				Ker = M.right_kernel()
				d = vector(QQ,M.ncols())
				for i in range(len(Ker.basis())):
					if w*Ker.basis_matrix()[i] > 0:
						d = Ker.basis_matrix()[i]
						break
					if w*Ker.basis_matrix()[i] < 0:
						d = -Ker.basis_matrix()[i]
						break
				print "LP is feasible, but unbounded "
				print "Check:"
				print M*(U.transpose()*x_ext + 10*d)<= v and w*(U.transpose()*x_ext + 10*d) > w*(U.transpose()*x_ext +d)
				return [U.transpose()*x_ext,d]
	else:
		# Decide whether LP feasible or not
		B = ind_rows(A,[])
		zero = vector(QQ,A.ncols())

		[x,l,B] = phase_II(A,B,b,zero)
		if len(x) > 0:
			x_ext = vector(QQ,M.ncols())
			d_ext = vector(QQ,M.ncols())
			for i in range(0,len(x)):
				x_ext[i] = x[i]
				d_ext[i] = d[i]
			print "LP is feasible, but unbounded"
			print "Check:"
			print M*(U.transpose()*x_ext + 10*U.transpose()*d_ext)<=v and w*(U.transpose()*x_ext + 10*U.transpose()*d_ext) >  w*(U.transpose()*x_ext + U.transpose()*d_ext)
			return [U.transpose()*x_ext,U.transpose()*d_ext]
		else:
			print "LP is infeasible"
			print "Check:"
			print l*M == 0 and l*v == -1
			return [x,l]


def randomTest (m, n):
    A = random_matrix (QQ, m, n)
    b = random_matrix (QQ, 1, m)[0]
    c = random_matrix (QQ, 1, n)[0]
#    print 'DEBUG: A:\n', A
#    print 'DEBUG: b: ', b
#    print 'DEBUG: c: ', c
    simplex(A,b,c)
    return