MecJan2025 update

.gitignore

@@ -1,2 +1,7 @@
 # OS or Editor folders
 .DS_Store
+.idea/.name
+.idea/mec_pypi.iml
+.idea/misc.xml
+.idea/modules.xml
+.idea/vcs.xml

mec/data.py

@@ -19,7 +19,7 @@ def load_stigler_data(nbi = 9, nbj = 77, verbose=False):
         print('\nDaily nutrient requirement:')
         print(allowance)
     return({'N_i_j':thedata.iloc[:nbj, 4:(4+nbi)].fillna(0).to_numpy().T,
-            'd_i':np.array(allowance)[0:nbi],
+            'b_i':np.array(allowance)[0:nbi],
             'c_j':np.ones(len(commodities))[0:nbj],
             'names_i': list(thedata.columns)[4:(4+nbi)],
             'names_j':commodities[0:nbj]}) 

mec/gt.py

@@ -36,7 +36,7 @@ def solve(self, verbose=0):
         return {'p_i': p_i, 'q_j': q_j, 'val': V_2}
 
     def simplex_solve(self, verbose=0):
-        tableau = Tableau(A_i_j = self.A_i_j, d_i = np.ones(self.nbi), c_j = np.ones(self.nbj),
+        tableau = Tableau(A_i_j = self.A_i_j, b_i = np.ones(self.nbi), c_j = np.ones(self.nbj),
                           decision_var_names_j = ['y_'+str(j) for j in range(self.nbj)])
         ystar, xstar, ystar_sum = tableau.simplex_solve()
         p_i = xstar / ystar_sum
@@ -133,7 +133,7 @@ def plot_cones(self):
 
     def create_tableau(self, display=False):
         tab = Tableau(A_i_j = -np.block([self.M_i_j, np.ones((self.nbi,1))]),
-                      d_i = self.q_i, c_j = None,
+                      b_i = self.q_i, c_j = None,
                       decision_var_names_j=self.z_names_i+['z_0'], slack_var_names_i=self.w_names_i)
         self.tableau = tab
         if display: tab.display()

mec/lp.py

@@ -30,13 +30,11 @@ def limited_tabulate(data, headers=None, tablefmt='grid', max_rows=18, max_cols=
 
 
 class LP():
-    def __init__(self, A_i_j, d_i, c_j = None, decision_var_names_j = None, slack_var_names_i = None):
-        if c_j is None:
-            c_j = np.zeros(A_i_j.shape[1])
+    def __init__(self, A_i_j, b_i, c_j = None, decision_var_names_j = None, slack_var_names_i = None):
         self.A_i_j = A_i_j
         self.nbi, self.nbj = A_i_j.shape
         self.nbk = self.nbi+self.nbj
-        self.d_i = d_i
+        self.b_i = b_i
         self.c_j = c_j
         if decision_var_names_j is None:
             decision_var_names_j = ['x_'+str(j) for j in range(self.nbj)]
@@ -51,7 +49,7 @@ def gurobi_solve(self,verbose=0):
             m.setParam('OutputFlag', 0)
         xg_j = m.addMVar(self.nbj)
         m.setObjective(xg_j@self.c_j,sense=grb.GRB.MAXIMIZE)
-        constr_i = m.addConstr(self.A_i_j @ xg_j <= self.d_i)
+        constr_i = m.addConstr(self.A_i_j @ xg_j <= self.b_i)
         m.optimize()
         return(xg_j.x,constr_i.pi,m.objVal)
 
@@ -60,21 +58,21 @@ def plot2d (self, the_path=[], legend=True):
         if len(self.c_j) != 2:
             print('The number of variables differs from two.')
             return()
-        x1max = min(di/self.A_i_j[i,0] for i, di in enumerate(self.d_i) if self.A_i_j[i,0] != 0 and di/self.A_i_j[i,0] >= 0)
-        x2max = min(di/self.A_i_j[i,1] for i, di in enumerate(self.d_i) if self.A_i_j[i,1] != 0 and di/self.A_i_j[i,1] >= 0)
+        x1max = min(di/self.A_i_j[i,0] for i, di in enumerate(self.b_i) if self.A_i_j[i,0] != 0 and di/self.A_i_j[i,0] >= 0)
+        x2max = min(di/self.A_i_j[i,1] for i, di in enumerate(self.b_i) if self.A_i_j[i,1] != 0 and di/self.A_i_j[i,1] >= 0)
         x1, x2 = np.meshgrid(np.linspace(-.2*x1max, 1.4*x1max, 400), np.linspace(-.2*x2max, 1.4*x2max, 400))
         feasible_region = (x1 >= 0) & (x2 >= 0)
-        for i, di in enumerate(self.d_i):
+        for i, di in enumerate(self.b_i):
             feasible_region = feasible_region & (self.A_i_j[i,0] * x1 + self.A_i_j[i,1] * x2 <= di)
         fig, ax = plt.subplots(figsize=(5, 5))
         plt.contourf(x1, x2, np.where(feasible_region, self.c_j[0]*x1 + self.c_j[1]*x2, np.nan), 50, alpha = 0.5, cmap='gray_r', levels=30)
-        for i, di in enumerate(self.d_i):
+        for i, di in enumerate(self.b_i):
             if self.A_i_j[i,1] != 0:
                 ax.plot(x1[0, :], di/self.A_i_j[i,1] - self.A_i_j[i,0]/self.A_i_j[i,1]*x1[0, :], label=self.slack_var_names_i[i]+' = 0')
             else:
                 ax.axvline(di/self.A_i_j[i,0], label=self.slack_var_names_i[i]+' = 0')
         if the_path:
-            ax.plot([a for (a,_) in the_path], [b for (_,b) in the_path], 'r--', label='Agorithm path')
+            ax.plot([a for (a,_) in the_path], [b for (_,b) in the_path], 'r--', label='Algorithm path')
             ax.scatter([a for (a,_) in the_path], [b for (_,b) in the_path], color='red')
         ax.set_xlim(-.2*x1max, 1.4*x1max), ax.set_ylim(-.2*x2max, 1.4*x2max)
         ax.set_xlabel(self.decision_var_names_j[0]), ax.set_ylabel(self.decision_var_names_j[1])
@@ -86,22 +84,28 @@ def plot2d (self, the_path=[], legend=True):
 
 
 class Dictionary(LP):
-    def __init__(self, A_i_j, d_i, c_j = None, slack_var_names_i=None,decision_var_names_j=None):
-        # s_i = d_i - A_i_j @ x_j
-        if d_i.min()<0:
+    def __init__(self, A_i_j, b_i, c_j = None, slack_var_names_i=None, decision_var_names_j=None):
+        # s_i = b_i - A_i_j @ x_j
+        if b_i.min() < 0:
             from warnings import warn
-            warn('The array d_i has negative entries; zero is not a feasible solution.')
-        LP.__init__(self, A_i_j, d_i, c_j, decision_var_names_j, slack_var_names_i)
+            warn('The array b_i has negative entries; zero is not a feasible solution.')
+        LP.__init__(self, A_i_j, b_i, c_j, decision_var_names_j, slack_var_names_i)
         self.nonbasic = [Symbol(x) for x in self.decision_var_names_j]
-        self.base = { Symbol('obj') : c_j @ self.nonbasic }
-        slack_exprs_i = d_i  - A_i_j @ self.nonbasic
+        if c_j is None: # for LCPs
+            self.base = {}
+        else: # for LPs
+            self.base = { Symbol('obj') : c_j @ self.nonbasic }
+        slack_exprs_i = b_i  - A_i_j @ self.nonbasic
         self.base.update({Symbol(name): slack_exprs_i[i] for (i,name) in enumerate(self.slack_var_names_i) })
 
     def variables(self):
-        return( list(self.base.keys())[1:] + self.nonbasic)
+        return( list(self.base.keys())[1:] + self.nonbasic )
 
     def display(self):
-        print('-------------------------- \nObjective and constraints:')
+        if self.c_j is None: # for LCPs
+            print('-------------------------- \nObjective and constraints:')
+        else: # for LPs
+            print('-------------------------- \nConstraints:')
         for var in self.base:
             print(var, '=', round_expr(self.base[var],2))
 
@@ -171,10 +175,10 @@ def dual_solution(self,verbose = 0):
 
 
     def simplex_loop(self,verbose = 0):
-        if self.d_i.min()<0:
+        if self.b_i.min()<0:
             from warnings import warn
-            warn('The array d_i has negative entries; zero is not a feasible solution.')
-        if verbose >2:
+            warn('The array b_i has negative entries; zero is not a feasible solution.')
+        if verbose>2:
             [x1,x2] = [Symbol(x) for x in self.decision_var_names_j]
             the_path = [self.primal_solution()]
         finished = False
@@ -185,14 +189,51 @@ def simplex_loop(self,verbose = 0):
         objVal = self.base[Symbol('obj')].subs([ (variable,0) for variable in self.nonbasic])
         if verbose>0:
             print('\nValue = ' + str(objVal))
-        if verbose >2:
+        if verbose>2:
             self.plot2d(the_path, legend=False)
         return (self.primal_solution(),self.dual_solution(),objVal)
 
+    def solution(self, verbose=0): # returns primal decision and slack variables
+        x_j = np.zeros(self.nbj)
+        s_i = np.zeros(self.nbi)
+        for j,var in enumerate(symbols(self.decision_var_names_j)):
+            if var in self.nonbasic:
+                x_j[j] = 0
+            else:
+                x_j[j] = float(self.base[var].subs([(variable,0) for variable in self.nonbasic]))
+            if verbose > 0: print(str(var) + ' = ' + str(x_j[j]))
+        for i,var in enumerate(symbols(self.slack_var_names_i)):
+            if var in self.nonbasic:
+                s_i[i] = 0
+            else:
+                s_i[i] = float(self.base[var].subs([(variable,0) for variable in self.nonbasic]))
+            if verbose > 0: print(str(var) + ' = ' + str(s_i[i]))
+        if self.c_j is not None: # for LPs
+            obj = float(self.base[Symbol('obj')].subs([(variable,0) for variable in self.nonbasic]))
+            if verbose > 0: print('obj' + ' = ' + str(obj))
+            return x_j, s_i, obj
+        else: # for LCPs
+            return x_j, s_i
+
+    def simplex_solve(self,verbose=0):
+        if self.b_i.min()<0:
+            from warnings import warn
+            warn('The array b_i has negative entries; zero is not a feasible solution.')
+        finished = False
+        x_sol, _, _ = self.solution()
+        path = [x_sol]
+        while not finished:
+            finished = self.step()
+            x_sol, _, _ = self.solution()
+            path.append(x_sol)
+        if verbose>0:
+            print('\nValue = '+str(example_dict2.base[Symbol('obj')].subs([ (variable,0) for variable in example_dict2.nonbasic])))
+        return self.solution()
+
 
 class Tableau(LP):
-    def __init__(self, A_i_j, d_i, c_j = None, slack_var_names_i = None, decision_var_names_j = None): # A_i_j @ x_j + s_i = d_i
-        LP.__init__(self, A_i_j, d_i, c_j, decision_var_names_j, slack_var_names_i)
+    def __init__(self, A_i_j, b_i, c_j = None, slack_var_names_i = None, decision_var_names_j = None): # A_i_j @ x_j + s_i = b_i
+        LP.__init__(self, A_i_j, b_i, c_j, decision_var_names_j, slack_var_names_i)
         self.nbi,self.nbj = A_i_j.shape
         self.nbk = self.nbi + self.nbj
         if c_j is None:
@@ -203,13 +244,14 @@ def __init__(self, A_i_j, d_i, c_j = None, slack_var_names_i = None, decision_va
             slack_var_names_i = ['s_'+str(i) for i in range(self.nbi)]
         self.names_all_variables =  self.slack_var_names_i + self.decision_var_names_j
         self.tableau = np.block( [[np.zeros((1,self.nbi)), c_j.reshape((1,-1)), 0],
-                                  [np.eye(self.nbi), A_i_j, d_i.reshape((-1,1))]] )
+                                  [np.eye(self.nbi), A_i_j, b_i.reshape((-1,1))]] )
         self.k_b = list(range(self.nbi)) # columns associated with basic variables
         self.i_b = list(range(1,1+self.nbi)) # rows associated with basic variables
 
     def display(self):
         tableau = []
-        tableau.append( ['Obj'] + list(self.tableau[0,:]) )
+        if c_j is not None:
+            tableau.append( ['Obj'] + list(self.tableau[0,:]) )
         for b in range(self.nbi):
             tableau.append([self.names_all_variables[self.k_b[b]]]+list(self.tableau[self.i_b[b],:]) )
         print(limited_tabulate(tableau, headers=[''] + self.names_all_variables + ['RHS'], tablefmt="grid"))
@@ -263,9 +305,9 @@ def simplex_step(self,verbose=0):
         return (kent is not None) # returns false  if optimal solution; true otherwise
 
     def simplex_solve(self,verbose=0):
-        if self.d_i.min()<0:
+        if self.b_i.min()<0:
             from warnings import warn
-            warn('The array d_i has negative entries; zero is not a feasible solution.')
+            warn('The array b_i has negative entries; zero is not a feasible solution.')
         while self.simplex_step(verbose):
             pass
         return self.solution()
@@ -321,12 +363,12 @@ def IP_loop(self, tol=1e-6, verbose=0):
         return True # finished
 
 
-def two_phase(A_i_j,d_i, verbose = False):
+def two_phase(A_i_j, b_i, verbose = False):
     nbi,nbj = A_i_j.shape
-    signs_i = np.minimum(2*np.sign(d_i)+1,1) # 1 if >=0, -1 else
-    d_i = signs_i*d_i
+    signs_i = np.minimum(2*np.sign(b_i)+1,1) # 1 if >=0, -1 else
+    b_i = signs_i * b_i
     A_i_j = signs_i[:,None] * A_i_j
-    the_tableau = Tableau(A_i_j, d_i, c_j = A_i_j.sum(axis= 0) )
+    the_tableau = Tableau(A_i_j, b_i, c_j = A_i_j.sum(axis= 0) )
     the_tableau.simplex_solve()
     if (min(the_tableau.k_b) >= nbi ):
         if verbose:

mec/nf.py

@@ -101,7 +101,7 @@ def __init__(self, nodesList, arcsList, c_a, q_z, active_basis=None, zero_node=0
         arcsNames = [str(x)+str(y) for (x,y) in self.arcsList]
 
         self.tableau = Tableau(A_i_j = np.asarray(np.linalg.solve(self.B(active_basis),self.N(active_basis))),
-                               d_i = self.q0_z,
+                               b_i = self.q0_z,
                                c_j = self.gain_a(active_basis)[self.nonbasis(active_basis)],
                                slack_var_names_i = [arcsNames[i] for i in active_basis],
                                decision_var_names_j = [arcsNames[i] for i in self.nonbasis(active_basis)])

mec/nf_newest.py

@@ -101,7 +101,7 @@ def __init__(self, nodesList, arcsList, c_a, q_z, active_basis=None, zero_node=0
         arcsNames = [str(x)+str(y) for (x,y) in self.arcsList]
 
         self.tableau = Tableau(A_i_j = np.asarray(np.linalg.solve(self.B(active_basis),self.N(active_basis))),
-                               d_i = self.q0_z,
+                               b_i = self.q0_z,
                                c_j = self.gain_a(active_basis)[self.nonbasis(active_basis)],
                                slack_var_names_i = [arcsNames[i] for i in active_basis],
                                decision_var_names_j = [arcsNames[i] for i in self.nonbasis(active_basis)])

mec/test.py

@@ -4,14 +4,14 @@
 def test_mec_lp_LP():
     from mec.lp import LP   
     data = mec.data.load_stigler_data(verbose = True)
-    dietLP = LP(data['N_i_j'].T,data['c_j'],data['d_i'])
+    dietLP = LP(data['N_i_j'].T,data['c_j'],data['b_i'])
     dietLP.gurobi_solve(verbose=0)
     return
 
 def test_mec_lp_Dictionary():
     from mec.lp import Dictionary
     dictionary_example = Dictionary(A_i_j = np.array([[2, 1], [1, 2]]),
-                                d_i = np.array([2,2]),
+                                b_i = np.array([2,2]),
                                 c_j = np.array([1,1]),
                                 slack_var_names_i = ['s_1', 's_2'],
                                 decision_var_names_j = ['x_1', 'x_2']