# =================== OPTIMISATION DISCRETE =================
# =                           EPFLove                       =
# ===========================================================

# Load simplex_corrige.py before executing the program.

# As we consider only matchings between a man and a woman, the corresponding graph is a complete bipartite graph.
# The weights on each edge in the graph is determined by the sum of the scores for the edge.
# We have seen in the lecture that the node-edge incidence matrix of a bipartite graph is totally unimodular.
# By the Hoffman-Kruskal theorem (p. 119 on the slides) we get that the linear relaxation of the maximum matching integer program
# is integral. This means that the simplex algorithm will find an integer optimum solution (Lemma p. 118).

def couplage(m,f,M,F):

	# Cree la matrice correspondant au systeme a resoudre
	# La taille pour les lignes correspondent a :
	# 		f lignes pour la condition somme sur delta(v) <= 1
	#		m lignes pour la condition somme sur delta(v) <= 1
	#		m*f lignes pour x_i >= 0
	# La taille des colonnes correspondent a :
	#		la correspondance respective avec chaque homme ou femme
	A = matrix(QQ, f+m+m*f, m*f)
	
	for i in range(f):
		for j in range(m):
			A[i,i*m + j] = 1
	
	for i in range(m):
		for j in range(f):
			A[i+f, i+m*j] = 1
			
					
	for i in range(m*f):
		A[m+f+i, i] = -1

	# Creation du vecteur c
	c = vector(QQ, m*f)
	
	for i in range(f):
		for j in range(m):
			c[i*m + j] = F[i, j] + M[j, i] 
	
	# Creation du vecteur b
	# Les f premiers composantes sont juste la condition somme sur delta(v) <= 1
	# Les m*f dernieres lignes sont pour la conditions -x_i <= 0
	b = vector(QQ, f+m+m*f)
	
	for i in range(f+m):
		b[i] = 1
	
	# On vient de creer un PL qui a des solutions entieres (cf nos explication du depart)
	# On peut alors appliquer notre algorithme du simplexe afin de maximiser tout ceci.
	[x,l] = simplex(A,b,c)
	
	print 'Le maximum de fric que peut se faire l\'agence est de : ' + str(c*x)
	print 'Les relations trouvees sont les suivantes : '
	for i in range(f):
		for j in range(m):
			if x[i*m + j] == 1:
				print '          Homme ' + str(j+1) + ' avec femme ' + str(i+1)

def couplage_stable(m,f,M,F):

	# Cree la matrice correspondant au systeme a resoudre
	# La taille pour les lignes correspondent a :
	# 		f lignes pour la condition somme sur delta(v) <= 1
	#		m lignes pour la condition somme sur delta(v) <= 1
	#		m*f lignes pour x_i >= 0
	#		m*f lignes pour les conditions de stabilite
	# La taille des colonnes correspondent a :
	#		la correspondance respective avec chaque homme ou femme
	A = matrix(QQ, f+m+m*f+m*f, m*f)
	
	for i in range(f):
		for j in range(m):
			A[i,i*m + j] = 1
	
	for i in range(m):
		for j in range(f):
			A[i+f, i+m*j] = 1		
					
	for i in range(m*f):
		A[m+f+i, i] = -1

	
	for i in range(m):
		for j in range(f):
			A[m*f + m + f + i+m*j, i+m*j] = -1
			for k in range(m):
				if F[j,k] > F[j,i]:
					A[m*f + m + f + i+m*j, k+m*j] = -1
				
			for l in range(f):
				if M[i,l] > M[i,j]:
					A[m*f + m + f + i+m*j, i+m*l] = -1

	# Creation du vecteur c
	c = vector(QQ, m*f)
	
	for i in range(f):
		for j in range(m):
			c[i*m + j] = F[i, j] + M[j, i] 
	
	# Creation du vecteur b
	# Les f premiers composantes sont juste la condition somme sur delta(v) <= 1
	# Les m*f dernieres lignes sont pour la conditions -x_i <= 0
	b = vector(QQ, f+m+m*f+m*f)
	
	for i in range(f+m):
		b[i] = 1
	
	for i in range(f*m):
		b[f+m+m*f+i] = -1

	# On vient de creer un PL qui a des solutions entieres (cf nos explication du depart)
	# On peut alors appliquer notre algorithme du simplexe afin de maximiser tout ceci.
	[x,l] = simplex(A,b,c)
	
	print 'Le maximum de fric que peut se asd faire l\'agence est de : ' + str(c*x)
	print 'Les relations trouvees sont les suivantes : '
	for i in range(f):
		for j in range(m):
			if x[i*m + j] == 1:
				print '          Homme ' + str(j+1) + ' avec femme ' + str(i+1)

# This variant uses the (MI)LP solver that SAGE has built in. It is in practice much faster than our simplex method.	
def couplage1(m,f,M,F):

	C = []
	p = MixedIntegerLinearProgram(maximization=True)
	x = p.new_variable(integer=False)

	profit = M + F.transpose()

	# set objective
	objective = 0
	for i in range(m):
		for j in range(f):
			objective += profit[i][j]*x[i*f+j]
	
	p.set_objective(objective)

	# set matching constraints:
	for i in range(m):
		constraint = 0
		for j in range(f):
			constraint += x[i*f+j]

		p.add_constraint(constraint <= 1)
	
	for j in range(f):
		constraint = 0
		for i in range(m):
			constraint += x[i*f+j]

		p.add_constraint(constraint <= 1)


	# stable matching constraints
		
	for i in range(m):
		for j in range(f):
			constraint = 0
			constraint += x[i*f +j]
			for k in range(m):
				if F[j,k] > F[j,i]:
					constraint += x[k*f + j]
				
			for l in range(f):
				if M[i,l] > M[i,j]:
					constraint += x[i*f +l]
			p.add_constraint(constraint >= 1)

	#p.show()

	objective = p.solve()

	for i, v in p.get_values(x).iteritems():
		#print 'x_%s = %s' % (i, v)
		if (v > 0.1):
			man = floor(i/f)
			woman = i%f
			C.append((man,woman))

	print objective
	return C

# Note that our condition of stability assumes that no two score values of a man (or a woman) are the same.
def random_agency(m,f):
	while True:
		M = matrix(QQ, m,f,[randint(0,100*m*f) for i in range(m*f)])
		if check_for_order(M):
			break
	while True:
		F = matrix(QQ, f,m,[randint(0,100*f*m) for i in range(m*f)])
		if check_for_order(F):
			break
	return [M,F]

def check_for_order(M):
	V = flatten (map (lambda x: list(x), M.rows ()))
	V.sort()
	for i in range(len(V) - 1):
		if V[i] == V[i+1]:
			return False
	return True