Skip to content

Commit

Permalink
v__0.0.8.9
Browse files Browse the repository at this point in the history 
added LCP class to the gt file
  • Loading branch information
@antoine-jacquet
antoine-jacquet committed Dec 14, 2023
1 parent b439133 commit 4226d62
Showing 1 changed file with 112 additions and 13 deletions.
125 changes: 112 additions & 13 deletions mec/gt.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,6 +56,107 @@ def minimax_CP(self,gap_threshold = 1e-5,max_iter = 10000):
return(p,q,gap,i)


class LCP: # z >= 0, w = M z + q >= 0, z.w = 0
def __init__(self,M_i_j,q_i):
self.M_i_j = M_i_j
self.nbi,_ = M_i_j.shape
self.q_i = q_i

def qp_solve(self, silent = True, verbose = 0):
qp = grb.Model()
if silent:
qp.Params.OutputFlag = 0
qp.Params.NonConvex = 2
zqp_i = qp.addMVar(shape = self.nbi)
wqp_i = qp.addMVar(shape = self.nbi)
qp.addConstr(self.M_i_j @ zqp_i - wqp_i == - self.q_i)
qp.setObjective(zqp_i @ wqp_i, sense = grb.GRB.MINIMIZE)
qp.optimize()
z_i = np.array(qp.getAttr('x'))[:self.nbi]
if verbose > 0: print(z_i)
if verbose > 1:
w_i = self.M_i_j @ z_i + self.q_i
print(w_i)
print(z_i @ w_i)
return(z_i)

def plot_cones(self):
if self.nbi != 2:
raise ValueError('Can\'t plot in 2D because the problem is of order different from 2.')
A, q = -self.M_i_j, self.q_i
I = np.eye(2)
fig = plt.figure(figsize=(5, 5))
ax = fig.add_subplot(1, 1, 1)
ax.spines['left'].set_position('zero'), ax.spines['bottom'].set_position('zero')
ax.spines['right'].set_color('none'), ax.spines['top'].set_color('none')
r = np.array([[.6,.7],[.8,.9]])
angles = np.array([[0, np.pi/2], [None, None]])
for j in range(2):
if any(A[:,j] != 0): angles[1,j] = (-1)**(A[1,j]<0) * np.arccos(A[0,j]/np.linalg.norm(A[:,j]))
for i in range(2):
for j in range(2):
angle_1 = angles[i,0]
angle_2 = angles[j,1]
if angle_1 != None and angle_2 != None:
if abs(angle_1 - angle_2) < np.pi:
theta = np.linspace(np.minimum(angle_1, angle_2), np.maximum(angle_1, angle_2), 100)
elif abs(angle_1 - angle_2) > np.pi:
theta = np.linspace(np.maximum(angle_1, angle_2), np.minimum(angle_1, angle_2) + 2*np.pi, 100)
ax.fill(np.append(0, r[i,j]*np.cos(theta)), np.append(0, r[i,j]*np.sin(theta)), 'b', alpha=0.1)
for i in range(2):
ax.quiver(0, 0, I[0,i], I[1,i], scale=1, scale_units='xy', color='b')
ax.text(I[0,i], I[1,i], r'$E_{{{}}}$'.format(i+1), fontsize=12, ha='left', va='bottom')
if any(A[:,i] != 0):
ax.quiver(0, 0, A[0,i]/np.linalg.norm(A[:,i]), A[1,i]/np.linalg.norm(A[:,i]), scale=1, scale_units='xy', color='b')
ax.text(A[0,i]/np.linalg.norm(A[:,i]), A[1,i]/np.linalg.norm(A[:,i]), r'$-M_{{\cdot{}}}$'.format(i+1), fontsize=12, ha='left', va='bottom')
ax.quiver(0, 0, q[0]/np.linalg.norm(q), q[1]/np.linalg.norm(q), angles='xy', scale_units='xy', scale=1, color='r')
plt.xlim(-1.2, 1.2), plt.ylim(-1.2, 1.2)
plt.show()

def lemke_solve(self,verbose=0,maxit=100):
counter = 0
zis = ['z_' + str(i+1) for i in range(self.nbi)]+['z_0']
wis = ['w_' + str(i+1) for i in range(self.nbi)]
if all(self.q_i >= 0):
print('==========')
print('Solution found: trivial LCP (q >= 0).')
zsol=np.zeros(self.nbi)
wsol=self.q_i
for i in range(self.nbi): print(zis[i] + ' = 0')
for i in range(self.nbi): print(wis[i] + ' = ' + str(wsol[i]))
return zsol,wsol
else:
keys = wis + zis[:-1]
labels = zis[:-1] + wis
complements = {Symbol(keys[i]): Symbol(labels[i]) for i in range(2*self.nbi)}
tab = Tableau(wis, zis, - np.block([self.M_i_j, np.ones((self.nbi,1))]), self.q_i)
entering_var = Symbol('z_0')
departing_var = Symbol(wis[np.argmin(self.q_i)]) # if argmin not unique, takes smallest index (Bland's rule)
while departing_var is not None and counter < maxit:
counter += 1
tab.pivot(entering_var,departing_var,verbose=verbose)
if departing_var == Symbol('z_0') or counter == maxit:
break
else:
entering_var = complements[departing_var]
departing_var = tab.determine_departing(entering_var)
print('==========')
if departing_var == Symbol('z_0'):
print('Solution found: z_0 departed basis.')
sol=tab.solution()
zsol=np.array([sol[Symbol(zis[i])] for i in range(self.nbi)])
wsol=np.array([sol[Symbol(wis[i])] for i in range(self.nbi)])
print('Complementarity: z.w = ' + str(zsol @ wsol))
for i in range(2*self.nbi): print(labels[i] + ' = ' + str(sol[Symbol(labels[i])]))
return zsol, wsol
elif counter == maxit:
print('Solution not found: maximum number of iterations (' + str(maxit) + ') reached.')
return None
else:
print('Solution not found: Ray termination.')
return None


class Bimatrix_game:
def __init__(self,A_i_j,B_i_j):
self.A_i_j = A_i_j
Expand All @@ -82,15 +183,18 @@ def mangasarian_stone_solve(self):

def lemke_howson_solve(self,verbose = 0):
from sympy import Symbol


A_i_j = self.A_i_j - np.min(self.A_i_j) + 1 # adding constant to ensure that matrices are positive
B_i_j = self.B_i_j - np.min(self.B_i_j) + 1

ris = ['r_' + str(i+1) for i in range(self.nbi)]
yjs = ['y_' + str(self.nbi+j+1) for j in range(self.nbj)]
sjs = ['s_' + str(self.nbi+j+1) for j in range(self.nbj)]
xis = ['x_' + str(i+1) for i in range(self.nbi)]
#tab2 = Tableau(ris, yjs, self.A_i_j, np.ones(self.nbi) )
tab2 = Dictionary( self.A_i_j, np.ones(self.nbi),np.zeros(self.nbj),ris, yjs )
#tab1 = Tableau(sjs, xis, self.B_i_j.T, np.ones(self.nbj) )
tab1 = Dictionary(self.B_i_j.T, np.ones(self.nbj), np.zeros(self.nbi), sjs, xis)
#tab2 = Tableau(ris, yjs, A_i_j, np.ones(self.nbi) )
tab2 = Dictionary(A_i_j, np.ones(self.nbi),np.zeros(self.nbj),ris, yjs )
#tab1 = Tableau(sjs, xis, B_i_j.T, np.ones(self.nbj) )
tab1 = Dictionary(B_i_j.T, np.ones(self.nbj), np.zeros(self.nbi), sjs, xis)
keys = ris+yjs+sjs+xis
labels = xis+sjs+yjs+ris
complements = {Symbol(keys[t]): Symbol(labels[t]) for t in range(len(keys))}
Expand All @@ -117,7 +221,7 @@ def lemke_howson_solve(self,verbose = 0):
val2 = 1 / x_i.sum()
p_i = x_i * val2
q_j = y_j * val1
sol_dict = {'val1':val1, 'val2':val2, 'p_i':p_i,'q_j':q_j}
sol_dict = {'val1': val1, 'val2': val2, 'p_i': p_i, 'q_j': q_j}
return(sol_dict)

import numpy as np
Expand Down Expand Up @@ -162,7 +266,6 @@ def is_standard_form(self):
cond_1 = (np.diag(self.Phi_z_a) == self.Phi_z_a.min(axis = 1) ).all()
cond_2 = ((self.Phi_z_a[:,:self.nbz] + np.diag([np.inf] * self.nbz)).min(axis=1) >= self.Phi_z_a[:,self.nbz:].max(axis=1)).all()
return (cond_1 & cond_2)


def p_z(self,basis=None):
if basis is None:
Expand All @@ -177,7 +280,6 @@ def musol_a(self,basis=None):
mu_a[list(basis)] = np.linalg.solve(B,self.q_z)
return mu_a


def is_feasible_basis(self,basis):
try:
if self.musol_a(list(basis) ).min()>=0:
Expand All @@ -193,8 +295,7 @@ def is_ordinal_basis(self,basis):
if blocking.any():
which = np.where(blocking)
return res, which



def determine_entering(self,a_departing):
self.nbstep += 1
pbefore_z = self.p_z(self.basis_Phi)
Expand All @@ -207,9 +308,7 @@ def determine_entering(self,a_departing):
a_entering = max([(c,self.Phi_z_a[zstar,c]) for c in eligible_columns], key = lambda x: x[1])[0]
self.basis_Phi.append(a_entering)
return a_entering




def step(self,a_entering ,verbose= 0):
a_departing = self.tableau_M.determine_departing(a_entering)
self.tableau_M.update(a_entering,a_departing)
Expand Down

0 comments on commit 4226d62

Please sign in to comment.