From 7eda9afb64fff240ca35c6b38be9645ebdc9b255 Mon Sep 17 00:00:00 2001 From: Rex-8 Date: Mon, 30 Mar 2026 20:57:29 +0530 Subject: [PATCH] refactor(cli): replace sys.argv parsing with click-based CLI --- src/lemke/bimatrix.py | 389 +++++++++++++------------- src/lemke/lemke.py | 614 +++++++++++++++++++----------------------- 2 files changed, 459 insertions(+), 544 deletions(-) diff --git a/src/lemke/bimatrix.py b/src/lemke/bimatrix.py index 01b329c..53b128f 100644 --- a/src/lemke/bimatrix.py +++ b/src/lemke/bimatrix.py @@ -1,23 +1,25 @@ # bimatrix class +import fractions +import random # random.seed import sys + +import click import numpy as np -import fractions -from . import utils -from . import columnprint -from . import lemke -from . import randomstart -import random # random.seed + +from . import columnprint, lemke, randomstart, utils + # for debugging -def printglobals(): +def printglobals(): globs = [x for x in globals().keys() if not "__" in x] for var in globs: value = str(globals()[var]) if not "<" in value: - print(" "+str(var)+"=", value) + print(" " + str(var) + "=", value) + -# file format: +# file format: # # m*n entries of A, separated by blanks / newlines # m*n entries of B, separated by blanks / newlines @@ -25,123 +27,27 @@ def printglobals(): # blank lines or lines starting with "#" are ignored # defaults -# MAXDIM = 2000 # largest allowed value for m and n; not used yet gamefilename = "game" gz0 = False -LHstring = "" # empty means LH not called +LHstring = "" # empty means LH not called seed = -1 -trace = -1 # negative: no tracing -# accuracy = DEFAULT_accuracy = 1000 +trace = -1 # negative: no tracing accuracy = 1000 -# amends defaults -def processArguments(): - global gamefilename,gz0,LHstring,seed,trace,accuracy - arglist = sys.argv[1:] - setLH = False - settrace = False - setaccuracy = False - setseed = False - setdecimals = False - showhelp = False - for s in arglist: - if s[0] == '-': - # discard optional argument-parameters - if (setLH or settrace or setseed): - setLH = settrace = setseed = False - if s=="-LH": - setLH = True - LHstring = "1-" - elif s=="-trace": - settrace = True - trace = 0 - elif s=="-seed": - setseed = True - seed = 0 - elif s=="-accuracy": - setaccuracy = True - elif s=="-decimals": - setdecimals = True - elif s=="-z0": - gz0 = True - elif s=="-?" or s=="-help": - showhelp = True - else: # any other "-" argument - showhelp = True - print("! unknown option: ", repr(s)) - else: # s not starting with "-" - if setLH: - setLH = False - LHstring = s - elif settrace: - settrace= False - trace = int(s) - elif setseed: - setseed= False - seed = int(s) - elif setdecimals: - setdecimals = False - utils.setdecimals(int(s)) - elif setaccuracy: - setaccuracy= False - accuracy = int(s) - else: - gamefilename = s - if (showhelp): - helpstring="""usage: bimatrix.py [options] -options: - here: """+repr(gamefilename)+ """, must not start with '-' - -LH [] : Lemke-Howson with missing labels, e.g. '1,3-5,7-' ('' = all) - -trace []: tracing procedure, no. of priors, 0 = centroid - -seed [] : random seed, default: None - -accuracy : accuracy prior, =denominator, here """+str(accuracy)+""" - -decimals : allowed payoff digits in input after decimal point, default 4 - -?, -help: show this help and exit""" - print(helpstring) - exit(0) - return - -# list generated from string s such as "1-3,10,4-7", all not -# larger than endrange (50 is arbitrary default) and at least 1 -def rangesplit(s,endrange=50): - result = [] - for part in s.split(','): - if part != "": - if '-' in part: - a, b = part.split('-') - a = int(a) - b = endrange if b=="" else int(b) - else: - a = int(part) - b = a - a = max(a,1) - b = min(b, endrange) # a > endrange means empty range - result.extend(range(a, b+1)) - return result - + # used for both A and B class payoffmatrix: - # create zero matrix of given dimensions - def __init__(self, m, n): - self.numrows = m - self.numcolumns = n - self.matrix = np.zeros( (m,n), dtype=fractions.Fraction) - self.negmatrix = np.zeros( (m,n), dtype=fractions.Fraction) - self.max = 0 - self.min = 0 - self.negshift = 0 - # create matrix from any numerical matrix def __init__(self, A): AA = np.array(A) - m,n = AA.shape + m, n = AA.shape self.numrows = m self.numcolumns = n - self.matrix = np.zeros( (m,n), dtype=fractions.Fraction) + self.matrix = np.zeros((m, n), dtype=fractions.Fraction) for i in range(m): for j in range(n): self.matrix[i][j] = utils.tofraction(AA[i][j]) - self.fullmaxmin() + self.fullmaxmin() def __str__(self): buf = columnprint.columnprint(self.numcolumns) @@ -153,58 +59,51 @@ def __str__(self): out += ", negshift= " + str(self.negshift) return out - def updatemaxmin(self, fromrow, fromcol): - m=self.numrows - n=self.numcolumns + def updatemaxmin(self, fromrow, fromcol): + m = self.numrows + n = self.numcolumns for i in range(fromrow, m): for j in range(fromcol, n): elt = self.matrix[i][j] self.max = max(self.max, elt) self.min = min(self.min, elt) - self.negshift = int(self.max)+1 - self.negmatrix = np.full((m,n),self.negshift,dtype=int)-self.matrix + self.negshift = int(self.max) + 1 + self.negmatrix = np.full((m, n), self.negshift, dtype=int) - self.matrix def fullmaxmin(self): self.max = self.matrix[0][0] self.min = self.matrix[0][0] - self.updatemaxmin(0,0) + self.updatemaxmin(0, 0) - # add full row, row must be of size n - def addrow(self, row): + def addrow(self, row): self.matrix = np.vstack([self.matrix, row]) self.numrows += 1 - self.updatemaxmin(self.numrows-1,0) + self.updatemaxmin(self.numrows - 1, 0) - # add full column, col must be of size m def addcolumn(self, col): self.matrix = np.column_stack([self.matrix, col]) self.numcolumns += 1 - self.updatemaxmin(0,self.numcolumns-1) - -class bimatrix: - # create A,B given m,n - def __init__(self, m, n): - self.A = payoffmatrix(m,n) - self.B = payoffmatrix(m,n) + self.updatemaxmin(0, self.numcolumns - 1) + - # create A,B from file +class bimatrix: def __init__(self, filename): lines = utils.stripcomments(filename) - # flatten into words words = utils.towords(lines) m = int(words[0]) n = int(words[1]) - needfracs = 2*m*n + needfracs = 2 * m * n if len(words) != needfracs + 2: - print("in bimatrix file "+repr(filename)+":") - print("m=",m,", n=",n,", need", - needfracs,"payoffs, got", len(words)-2) + print("in bimatrix file " + repr(filename) + ":") + print( + "m=", m, ", n=", n, ", need", needfracs, "payoffs, got", len(words) - 2 + ) exit(1) k = 2 - C = utils.tomatrix(m, n, words, k) + C = utils.tomatrix(m, n, words, k) self.A = payoffmatrix(C) - k+= m*n - C = utils.tomatrix(m, n, words, k) + k += m * n + C = utils.tomatrix(m, n, words, k) self.B = payoffmatrix(C) def __str__(self): @@ -217,60 +116,57 @@ def __str__(self): def createLCP(self): m = self.A.numrows n = self.A.numcolumns - lcpdim = m+n+2 + lcpdim = m + n + 2 lcp = lemke.lcp(lcpdim) - lcp.q[lcpdim-2] = -1 - lcp.q[lcpdim-1] = -1 + lcp.q[lcpdim - 2] = -1 + lcp.q[lcpdim - 1] = -1 for i in range(m): - lcp.M[lcpdim-2][i] = 1 - lcp.M[i][lcpdim-2] = -1 - for j in range(m,m+n): - lcp.M[lcpdim-1][j] = 1 - lcp.M[j][lcpdim-1] = -1 + lcp.M[lcpdim - 2][i] = 1 + lcp.M[i][lcpdim - 2] = -1 + for j in range(m, m + n): + lcp.M[lcpdim - 1][j] = 1 + lcp.M[j][lcpdim - 1] = -1 for i in range(m): for j in range(n): - lcp.M[i][j+m] = self.A.negmatrix[i][j] + lcp.M[i][j + m] = self.A.negmatrix[i][j] for j in range(n): for i in range(m): - lcp.M[j+m][i] = self.B.negmatrix[i][j] - # d for now + lcp.M[j + m][i] = self.B.negmatrix[i][j] for i in range(lcpdim): - lcp.d[i]=1 + lcp.d[i] = 1 return lcp - def runLH(self, droppedlabel): lcp = self.createLCP() - lcp.d[droppedlabel-1] = 0 # subsidize this label + lcp.d[droppedlabel - 1] = 0 tabl = lemke.tableau(lcp) - # tabl.runlemke(verbose=True, lexstats=True, z0=gz0) tabl.runlemke(silent=True) return tuple(getequil(tabl)) - + def LH(self, LHstring): if LHstring == "": return m = self.A.numrows n = self.A.numcolumns - lhset = {} # dict of equilibria and list by which label found - labels = rangesplit(LHstring, m+n) + lhset = {} + labels = rangesplit(LHstring, m + n) for k in labels: eq = self.runLH(k) if eq in lhset: lhset[eq].append(k) else: - print ("label",k,"found eq", str_eq(eq,m,n)) - lhset[eq] = [k] - print ("-------- equilibria found: --------") + print("label", k, "found eq", str_eq(eq, m, n)) + lhset[eq] = [k] + print("-------- equilibria found: --------") for eq in lhset: - print (str_eq(eq,m,n),"found by labels", str(lhset[eq])) + print(str_eq(eq, m, n), "found by labels", str(lhset[eq])) return lhset def runtrace(self, xprior, yprior): lcp = self.createLCP() Ay = self.A.negmatrix @ yprior - xB = xprior @ self.B.negmatrix - lcp.d = np.hstack((Ay,xB,[1,1])) + xB = xprior @ self.B.negmatrix + lcp.d = np.hstack((Ay, xB, [1, 1])) tabl = lemke.tableau(lcp) tabl.runlemke(silent=True) return tuple(getequil(tabl)) @@ -280,85 +176,170 @@ def tracing(self, trace): return m = self.A.numrows n = self.A.numcolumns - trset = {} # dict of equilibria, how often found + trset = {} if trace == 0: xprior = uniform(m) yprior = uniform(n) eq = self.runtrace(xprior, yprior) - trset[eq]=1 - trace = 1 # for percentage + trset[eq] = 1 + trace = 1 else: for k in range(trace): - if seed >=0: - random.seed(10*trace*seed+k) + if seed >= 0: + random.seed(10 * trace * seed + k) x = randomstart.randInSimplex(m) xprior = randomstart.roundArray(x, accuracy) y = randomstart.randInSimplex(n) yprior = randomstart.roundArray(y, accuracy) - # print (f"{k=} {xprior=} {yprior=}") eq = self.runtrace(xprior, yprior) if eq in trset: trset[eq] += 1 else: - print ("found eq", str_eq(eq,m,n), "index", - self.eqindex(eq,m,n)) - trset[eq] = 1 - print ("-------- statistics of equilibria found: --------") + print("found eq", str_eq(eq, m, n), "index", self.eqindex(eq, m, n)) + trset[eq] = 1 + print("-------- statistics of equilibria found: --------") for eq in trset: - print (trset[eq],"times found ",str_eq(eq,m,n)) - print(trace,"total priors,",len(trset),"equilibria found") + print(trset[eq], "times found ", str_eq(eq, m, n)) + print(trace, "total priors,", len(trset), "equilibria found") - def eqindex(self,eq,m,n): - rowset,colset = supports(eq,m,n) - k,l = len(rowset),len(colset) - if k!=l: + def eqindex(self, eq, m, n): + rowset, colset = supports(eq, m, n) + k, l = len(rowset), len(colset) + if k != l: return 0 A1 = submatrix(self.A.negmatrix, rowset, colset) DA = np.linalg.det(A1) B1 = submatrix(self.B.negmatrix, rowset, colset) DB = np.linalg.det(B1) - sign = 2*(k%2) - 1 # -1 if even, 1 if odd - if DA*DB == 0: + sign = 2 * (k % 2) - 1 + if DA * DB == 0: return 0 - if DA*DB > 0: + if DA * DB > 0: return sign - return -sign - + return -sign + + def uniform(n): - return np.array([ fractions.Fraction(1,n) for j in range(n)]) + return np.array([fractions.Fraction(1, n) for j in range(n)]) + def getequil(tabl): tabl.createsol() - return tabl.solution[1:tabl.n-1] - -def str_eq(eq,m,n): - x = "("+",".join([str(x) for x in eq[0:m]])+")" - y = "("+",".join([str(x) for x in eq[m:m+n]])+")" - rowset,colset = supports(eq,m,n) - return x+","+y+"\n supports: "+str(rowset)+str(colset) - -def supports(eq,m,n): - rowset = [i for i in range(m) if eq[i]!= 0] - colset = [j for j in range(n) if eq[m+j]!= 0] - return rowset,colset - -def submatrix(A,rowset,colset): - k,l = len(rowset),len(colset) - B = np.zeros((k,l)) - for i in range (k): + return tabl.solution[1 : tabl.n - 1] + + +def str_eq(eq, m, n): + x = "(" + ",".join([str(x) for x in eq[0:m]]) + ")" + y = "(" + ",".join([str(x) for x in eq[m : m + n]]) + ")" + rowset, colset = supports(eq, m, n) + return x + "," + y + "\n supports: " + str(rowset) + str(colset) + + +def supports(eq, m, n): + rowset = [i for i in range(m) if eq[i] != 0] + colset = [j for j in range(n) if eq[m + j] != 0] + return rowset, colset + + +def submatrix(A, rowset, colset): + k, l = len(rowset), len(colset) + B = np.zeros((k, l)) + for i in range(k): for j in range(l): B[i][j] = A[rowset[i]][colset[j]] return B -def main(): - processArguments() +def rangesplit(s, endrange=50): + result = [] + for part in s.split(","): + if part != "": + if "-" in part: + a, b = part.split("-") + a = int(a) + b = endrange if b == "" else int(b) + else: + a = int(part) + b = a + a = max(a, 1) + b = min(b, endrange) + result.extend(range(a, b + 1)) + return result + + +# Sentinel: distinguishes "option not passed" from "option passed without value" +_UNSET = object() + + +@click.command() +@click.argument("gamefilename", default="game") +@click.option( + "-LH", + "lh_range", + default=None, + metavar="RANGE", + help="Lemke-Howson with missing labels, e.g. '1,3-5,7-' (omit value = all labels).", +) +@click.option( + "-trace", + "trace", + default=None, + type=int, + metavar="NUM", + help="Tracing procedure; NUM = number of priors (0 = centroid).", +) +@click.option( + "-seed", "seed_val", default=None, type=int, metavar="NUM", help="Random seed." +) +@click.option( + "-accuracy", + "accuracy_val", + default=1000, + type=int, + metavar="N", + help="Accuracy prior denominator (default: 1000).", +) +@click.option( + "-decimals", + "decimals_val", + default=None, + type=int, + metavar="D", + help="Allowed payoff digits after decimal point (default: 4).", +) +@click.option( + "-z0", "gz0_flag", is_flag=True, default=False, help="Use z0 entering variable." +) +def main(gamefilename, lh_range, trace, seed_val, accuracy_val, decimals_val, gz0_flag): + global gz0, LHstring, seed, accuracy + + # -LH with no argument in original code defaulted to "1-" (all labels). + # With Click, omitting -LH entirely means lh_range is None (LH not called). + # Passing -LH with an explicit range string works as before. + # To replicate "-LH with no argument" pass -LH "1-". + LHstring = lh_range if lh_range is not None else "" + + # -trace with no argument in original code defaulted to 0 (centroid). + # With Click, omitting -trace means trace stays -1 (tracing not called). + # Passing -trace 0 explicitly replicates the original "-trace" with no arg. + trace_val = trace if trace is not None else -1 + + # -seed with no argument in original code defaulted to 0. + # With Click, omitting -seed means seed stays -1 (no fixed seed). + seed = seed_val if seed_val is not None else -1 + + accuracy = accuracy_val + gz0 = gz0_flag + + if decimals_val is not None: + utils.setdecimals(decimals_val) + printglobals() G = bimatrix(gamefilename) print(G) - eqset = G.LH(LHstring) - eqset = G.tracing(trace) + G.LH(LHstring) + G.tracing(trace_val) if __name__ == "__main__": diff --git a/src/lemke/lemke.py b/src/lemke/lemke.py index 7c10e08..4431996 100644 --- a/src/lemke/lemke.py +++ b/src/lemke/lemke.py @@ -1,126 +1,98 @@ # LCP solver +import fractions +import math # gcd import sys -import fractions -import math # gcd -from . import columnprint -from . import utils +import click + +from . import columnprint, utils # global defaults -lcpfilename="lcp" -outfile=lcpfilename+".out" +lcpfilename = "lcp" +outfile = lcpfilename + ".out" filehandle = sys.stdout -verbose=False -silent=False -z0=False - -# process command-line arguments -def processArguments(): - global lcpfilename,outfile,filehandle,verbose,silent,z0 - helpstring="""usage: lemke.py [options] -options: -v, -verbose : printout intermediate tableaus - -s, -silent : send output to .out - -z0 : show value of z0 at each step - -?, -help: show this help - (default: "lcp", must not start with "-") - Example: [python] lemke.py -v 2lcp""" - arglist = sys.argv[1:] - showhelp = False - for s in arglist: - if s=="-v" or s=="-verbose": - verbose = True - elif s=="-s" or s=="-silent": - silent = True - elif s=="-z0": - z0 = True - elif s[0]=="-" : - showhelp = True - else: - lcpfilename = s - outfile = s+".out" - if (showhelp): - printout(helpstring) - exit(0) - return +verbose = False +silent = False +z0 = False + -def printout(*s): - print(*s, file=filehandle) - # LCP data M,q,d -class lcp: +class lcp: # create LCP either with given n or from file - def __init__(self, arg): - if isinstance(arg, int): # arg is an integer + def __init__(self, arg): + if isinstance(arg, int): # arg is an integer n = self.n = arg - # self.M = np.zeros( (n,n), dtype=fractions.Fraction) - # self.q = np.zeros( (n), dtype=fractions.Fraction) - # self.d = np.zeros( (n), dtype=fractions.Fraction) - self.M = [[]]*n + self.M = [[]] * n for i in range(n): - self.M[i]=[0]*n - self.q = [0]*n - self.d = [0]*n - else: # assume arg is a string = name of lcp file + self.M[i] = [0] * n + self.q = [0] * n + self.d = [0] * n + else: # assume arg is a string = name of lcp file # create LCP from file filename = arg lines = utils.stripcomments(filename) # flatten into words words = utils.towords(lines) - if words[0]!="n=": - printout("lcp file",repr(filename), - "must start with 'n=' lcpdim, e.g. 'n= 5', not", - repr(words[0])) + if words[0] != "n=": + printout( + "lcp file", + repr(filename), + "must start with 'n=' lcpdim, e.g. 'n= 5', not", + repr(words[0]), + ) exit(1) n = int(words[1]) self.n = n - # self.M = np.zeros( (n,n), dtype=fractions.Fraction) - # self.d = np.zeros( (n), dtype=fractions.Fraction) - # self.q = np.zeros( (n), dtype=fractions.Fraction) - self.M = [[]]*n + self.M = [[]] * n for i in range(n): - self.M[i]=[0]*n - self.q = [0]*n - self.d = [0]*n - needfracs = n*n + 2*n + self.M[i] = [0] * n + self.q = [0] * n + self.d = [0] * n + needfracs = n * n + 2 * n if len(words) != needfracs + 5: - # printout("in lcp file '",filename,"':") - printout("in lcp file "+repr(filename)+":") - printout("n=",n,", need keywords 'M=' 'q=' 'd=' and n*n + n + n =", - needfracs,"fractions, got", len(words)-5) + printout("in lcp file " + repr(filename) + ":") + printout( + "n=", + n, + ", need keywords 'M=' 'q=' 'd=' and n*n + n + n =", + needfracs, + "fractions, got", + len(words) - 5, + ) exit(1) - k = 2 # index in words + k = 2 # index in words while k < len(words): - if words[k]=="M=": - k+=1 - self.M = utils.tomatrix(n,n,words,k) - k+= n*n - elif words[k]=="q=": - k+=1 - self.q = utils.tovector(n,words,k) - k+=n - elif words[k]=="d=": - k+=1 - self.d = utils.tovector(n,words,k) - k+=n - else: - printout("in lcp file "+repr(filename)+":") - printout("expected one of 'M=' 'q=' 'd=', got",repr(words[k])) + if words[k] == "M=": + k += 1 + self.M = utils.tomatrix(n, n, words, k) + k += n * n + elif words[k] == "q=": + k += 1 + self.q = utils.tovector(n, words, k) + k += n + elif words[k] == "d=": + k += 1 + self.d = utils.tovector(n, words, k) + k += n + else: + printout("in lcp file " + repr(filename) + ":") + printout("expected one of 'M=' 'q=' 'd=', got", repr(words[k])) exit(1) - return - + return + def __str__(self): - n=self.n - M=self.M - q=self.q - d=self.d + n = self.n + M = self.M + q = self.q + d = self.d m = columnprint.columnprint(n) m.makeLeft(0) m.sprint("M=") m.newline() for i in range(n): - for j in range(n): - m.sprint(str(M[i][j])) + for j in range(n): + m.sprint(str(M[i][j])) m.sprint("q=") m.newline() for i in range(n): @@ -129,317 +101,271 @@ def __str__(self): m.newline() for i in range(n): m.sprint(str(d[i])) - # printout("M[0][0]", type(M[0][0])) - return "n= "+str(n)+"\n"+str(m) + return "n= " + str(n) + "\n" + str(m) + ####### end of class lcp - + + class tableau: # filling the tableau from the LCP instance Mqd - def __init__(self, Mqd): + def __init__(self, Mqd): self.n = Mqd.n n = self.n - self.scalefactor = [0]*(n+2) # 0 for z0, n+1 for RHS - # A = tableau, long integer entries - # self.A = np.zeros( (n,n+2), dtype=object) - self.A = [[]]*n + self.scalefactor = [0] * (n + 2) # 0 for z0, n+1 for RHS + self.A = [[]] * n for i in range(n): - self.A[i]=[0]*(n+2) + self.A[i] = [0] * (n + 2) self.determinant = 1 - self.lextested = [0]*(n+1) - self.lexcomparisons = [0]*(n+1) + self.lextested = [0] * (n + 1) + self.lexcomparisons = [0] * (n + 1) self.pivotcount = 0 - self.solution = [fractions.Fraction(0)]*(2*n+1) # all vars - # variable encodings: VARS = 0..2n = Z(0) .. Z(n) W(1) .. W(n) - # tableau columns: RHS n+1 - # bascobas[v] in 0..n-1: basic, bascobas[v] = tableau row - # bascobas[v] in n..2n: cobasic, bascobas[v]-n = tableau col - self.bascobas = [0]*(2*n+1) - # whichvar inverse of bascobas, shows which basic/cobasic vars - self.whichvar = [0]*(2*n+1) - for i in range(n+1): # variables Z(i) all cobasic - self.bascobas[i] = n+i - self.whichvar[n+i] = i - for i in range(n): # variables W(i+1) all basic - self.bascobas[n+1+i] = i - self.whichvar[i] = n+1+i + self.solution = [fractions.Fraction(0)] * (2 * n + 1) # all vars + self.bascobas = [0] * (2 * n + 1) + self.whichvar = [0] * (2 * n + 1) + for i in range(n + 1): # variables Z(i) all cobasic + self.bascobas[i] = n + i + self.whichvar[n + i] = i + for i in range(n): # variables W(i+1) all basic + self.bascobas[n + 1 + i] = i + self.whichvar[i] = n + 1 + i # determine scale factors, lcm of denominators - for j in range(n+2): + for j in range(n + 2): factor = 1 for i in range(n): - if j==0: + if j == 0: den = Mqd.d[i].denominator - elif j==n+1: # RHS + elif j == n + 1: # RHS den = Mqd.q[i].denominator else: - den = Mqd.M[i][j-1].denominator + den = Mqd.M[i][j - 1].denominator # least common multiple - factor *= den // math.gcd(factor,den) + factor *= den // math.gcd(factor, den) self.scalefactor[j] = factor # fill in column j of A for i in range(n): - if j==0: + if j == 0: den = Mqd.d[i].denominator num = Mqd.d[i].numerator - elif j==n+1: # RHS + elif j == n + 1: # RHS den = Mqd.q[i].denominator num = Mqd.q[i].numerator else: - den = Mqd.M[i][j-1].denominator - num = Mqd.M[i][j-1].numerator - self.A[i][j] = (factor//den) * num + den = Mqd.M[i][j - 1].denominator + num = Mqd.M[i][j - 1].numerator + self.A[i][j] = (factor // den) * num self.determinant = -1 return - + def __str__(self): - out = "Determinant: "+str(self.determinant) + out = "Determinant: " + str(self.determinant) n = self.n - tabl = columnprint.columnprint(n+3) + tabl = columnprint.columnprint(n + 3) tabl.makeLeft(0) - tabl.sprint("var") # headers - for j in range(n+1): - tabl.sprint(self.vartoa(self.whichvar[j+n])) + tabl.sprint("var") # headers + for j in range(n + 1): + tabl.sprint(self.vartoa(self.whichvar[j + n])) tabl.sprint("RHS") - tabl.sprint("scfa") # scale factors - for j in range(n+2): - if j == n+1: # RHS - tabl.sprint(str(self.scalefactor[n+1])) - elif self.whichvar[j+n] > n: # col j is some W + tabl.sprint("scfa") # scale factors + for j in range(n + 2): + if j == n + 1: # RHS + tabl.sprint(str(self.scalefactor[n + 1])) + elif self.whichvar[j + n] > n: # col j is some W tabl.sprint("1") else: - tabl.sprint(str(self.scalefactor[self.whichvar[j+n]])) - tabl.newline() # blank line + tabl.sprint(str(self.scalefactor[self.whichvar[j + n]])) + tabl.newline() # blank line for i in range(n): tabl.sprint(self.vartoa(self.whichvar[i])) - for j in range(n+2): + for j in range(n + 2): s = str(self.A[i][j]) - if s == "0" : - s = "." # replace 0 by dot + if s == "0": + s = "." # replace 0 by dot tabl.sprint(s) - out += "\n"+ str(tabl) - out += "\n"+ "-----------------end of tableau-----------------" + out += "\n" + str(tabl) + out += "\n" + "-----------------end of tableau-----------------" return out - - def vartoa(self, v): # variable as as string w1..wn or z0..zn - if (v > self.n): - return "w"+str(v-self.n) + + def vartoa(self, v): # variable as as string w1..wn or z0..zn + if v > self.n: + return "w" + str(v - self.n) else: - return "z"+str(v) - - def createsol(self): # get solution from current tableau + return "z" + str(v) + + def createsol(self): # get solution from current tableau n = self.n - for i in range(2*n+1): + for i in range(2 * n + 1): row = self.bascobas[i] - if row < n: # i is a basic variable - num = self.A[row][n+1] - # value of Z(i): scfa[Z(i)]*rhs[row] / (scfa[RHS]*det) - # value of W(i-n): rhs[row] / (scfa[RHS]*det) - if i <= n: # computing Z(i) - num *= self.scalefactor[i] - self.solution[i] = fractions.Fraction(num, - self.determinant*self.scalefactor[n+1]) - else: # i is nonbasic - self.solution[i]=fractions.Fraction(0) - - def outsol(self): # string giving solution, after createsol() - # printout in columns to check complementarity + if row < n: # i is a basic variable + num = self.A[row][n + 1] + if i <= n: # computing Z(i) + num *= self.scalefactor[i] + self.solution[i] = fractions.Fraction( + num, self.determinant * self.scalefactor[n + 1] + ) + else: # i is nonbasic + self.solution[i] = fractions.Fraction(0) + + def outsol(self): # string giving solution, after createsol() n = self.n - sol = columnprint.columnprint(n+2) + sol = columnprint.columnprint(n + 2) sol.sprint("basis=") - for i in range(n+1): - if (self.bascobas[i]0 and self.bascobas[n+i] 0 and self.bascobas[n + i] < n: # W(i) is a basic variable + s = self.vartoa(n + i) else: s = " " sol.sprint(s) sol.sprint("z=") - for i in range(2*n+1): + for i in range(2 * n + 1): sol.sprint(str(self.solution[i])) - if i == n: # new line since printouting slack vars w next - sol.sprint ("w=") - sol.sprint ("") # no W(0) - return str(sol) - - def assertbasic(self, v, info): # assert that v is basic - if (self.bascobas[v] >= self.n): - printout (info, "Cobasic variable", self.vartoa(v), - "should be basic") + if i == n: # new line since printouting slack vars w next + sol.sprint("w=") + sol.sprint("") # no W(0) + return str(sol) + + def assertbasic(self, v, info): # assert that v is basic + if self.bascobas[v] >= self.n: + printout(info, "Cobasic variable", self.vartoa(v), "should be basic") exit(1) return - def assertcobasic(self, v, info): # assert that v is cobasic - if (self.bascobas[v] < self.n): - printout (info, "Cobasic variable", self.vartoa(v), - "should be cobasic") + def assertcobasic(self, v, info): # assert that v is cobasic + if self.bascobas[v] < self.n: + printout(info, "Cobasic variable", self.vartoa(v), "should be cobasic") exit(1) return - def docupivot(self, leave, enter): # leave, enter in VARS - self.assertbasic (leave, "docupivot") - self.assertcobasic (enter, "docupivot") + def docupivot(self, leave, enter): # leave, enter in VARS + self.assertbasic(leave, "docupivot") + self.assertcobasic(enter, "docupivot") s = "leaving: " + self.vartoa(leave).ljust(5) s += "entering: " + self.vartoa(enter) - printout (s) - return + printout(s) + return def raytermination(self, enter): - printout("Ray termination when trying to enter",self.vartoa(enter)) - printout (self) + printout("Ray termination when trying to enter", self.vartoa(enter)) + printout(self) printout("Current basis not an LCP solution:") self.createsol() printout(self.outsol()) exit(1) - def testtablvars(self): # msg only if error, continue + def testtablvars(self): # msg only if error, continue n = self.n - for i in range(2*n+1): - if self.bascobas[self.whichvar[i]] != i : - # injective suffices - for j in range(2*n+1): - if j==i: - printout ("First problem for j=",j,":") - # printout (f"{j=} {self.bascobas[j]=} {self.whichvar[j]=}") - printout (f"j={j} self.bascobas[j]={self.bascobas[j]} self.whichvar[j]={self.whichvar[j]}") + for i in range(2 * n + 1): + if self.bascobas[self.whichvar[i]] != i: + for j in range(2 * n + 1): + if j == i: + printout("First problem for j=", j, ":") + printout( + f"j={j} self.bascobas[j]={self.bascobas[j]} self.whichvar[j]={self.whichvar[j]}" + ) break - return - - def complement(self, v): # Z(i),W(i) are complements + return + + def complement(self, v): # Z(i),W(i) are complements n = self.n if v == 0: - printout ("Attempt to find complement of z0") + printout("Attempt to find complement of z0") exit(1) if v > n: - return v-n + return v - n else: - return v+n + return v + n - # output statistics of minimum ratio test def outstatistics(self): n = self.n lext = self.lextested - stats = columnprint.columnprint(n+2) + stats = columnprint.columnprint(n + 2) stats.makeLeft(0) stats.sprint("lex-column") - for i in range(n+1): + for i in range(n + 1): stats.iprint(i) stats.sprint("times tested") - for i in range(n+1): + for i in range(n + 1): stats.iprint(lext[i]) - if lext[0]>0: # otherwise never a degeneracy + if lext[0] > 0: # otherwise never a degeneracy stats.sprint("% of pivots") - for i in range(0,n+1): - stats.iprint(round(lext[i]*100/self.pivotcount)) + for i in range(0, n + 1): + stats.iprint(round(lext[i] * 100 / self.pivotcount)) stats.sprint("avg comparisons") - for i in range(n+1): - if lext[i]>0: - x = round(self.lexcomparisons[i]*10/lext[0]) - stats.sprint(str(x/10.0)) + for i in range(n + 1): + if lext[i] > 0: + x = round(self.lexcomparisons[i] * 10 / lext[0]) + stats.sprint(str(x / 10.0)) else: stats.sprint("-") - printout(stats) - - # returns leave,z0leave - # leave = leaving variable in VARS, given by lexmin row, - # when enter in VARS is entering variable - # only positive entries of entering column tested. - # Boolean z0leave indicates that z0 can leave the - # basis, but the lex-minratio test is performed fully, - # so leave might not be the index of z0 + printout(stats) + def lexminvar(self, enter): n = self.n - A = self.A + A = self.A self.assertcobasic(enter, "Lexminvar") - col = self.bascobas[enter]-n # entering tableau column + col = self.bascobas[enter] - n # entering tableau column leavecand = [] # candidates(=rows) for leaving var - for i in range(n): # start with positives in entering col + for i in range(n): # start with positives in entering col if A[i][col] > 0: leavecand.append(i) if leavecand == []: self.raytermination(enter) - if len(leavecand)==1: # single positive entering value - z0leave = self.bascobas[0] == leavecand[0] - ## omitted from statistics: only one possible row - ## means no min-ratio test needed for leaving variable - ## self.lextested[0] += 1 - ## self.lexcomparisons[0] += 1 - - # as long as there is more than one leaving candidate, - # perform a minimum ratio test for the columns - # j in RHS,W(1)..W(n) in the tableau. - # That test has an easy known result if the test - # column is basic, or equal to the entering variable. - j = 0 # going through j = 0..n - while len(leavecand)>1: - if j > n: # impossible, perturbed RHS should have full rank + if len(leavecand) == 1: + z0leave = self.bascobas[0] == leavecand[0] + j = 0 # going through j = 0..n + while len(leavecand) > 1: + if j > n: printout("lex-minratio test failed") exit(1) self.lextested[j] += 1 self.lexcomparisons[j] += len(leavecand) - if j==0: - testcol = n+1 # RHS + if j == 0: + testcol = n + 1 # RHS else: - testcol = self.bascobas[n+j]-n # tabl col of W(j) - if testcol != col: # otherwise nothing changed - if testcol >= 0: - # not a basic testcolumn: perform minimum ratio tests - newcand = [ leavecand[0] ] - # newcand contains the new candidates - for i in range(1,len(leavecand)): - # investigate remaining candidates - # compare ratios via products + testcol = self.bascobas[n + j] - n # tabl col of W(j) + if testcol != col: + if testcol >= 0: + newcand = [leavecand[0]] + for i in range(1, len(leavecand)): tmp1 = A[newcand[0]][testcol] * A[leavecand[i]][col] tmp2 = A[leavecand[i]][testcol] * A[newcand[0]][col] - # sgn = np.sign(tmp1 - tmp2) - # if sgn==0: - if tmp1 == tmp2 : # new ratio is the same as before + if tmp1 == tmp2: newcand.append(leavecand[i]) - elif tmp1 > tmp2: # new smaller ratio detected: reset - newcand = [ leavecand[i] ] - # else : unchanged candidates + elif tmp1 > tmp2: + newcand = [leavecand[i]] leavecand = newcand - else: # testcol < 0: W(j) basic, eliminate its row - # from leavecand if in there, since testcol is - # the jth unit column (ratio too big) - wj = self.bascobas[j+n] + else: + wj = self.bascobas[j + n] if wj in leavecand: leavecand.remove(wj) - # end of if testcol != col - # check if z0 among the first-col leaving candidates if j == 0: z0leave = self.bascobas[0] in leavecand - j += 1 # end while - assert (len(leavecand)==1) + j += 1 + assert len(leavecand) == 1 return self.whichvar[leavecand[0]], z0leave - # end of lexminvar(enter) - # negate tableau column col def negcol(self, col): for i in range(self.n): - self.A[i][col] = -self.A[i][col] - - # negate tableau row. Used in pivot() + self.A[i][col] = -self.A[i][col] + def negrow(self, row): - for j in range(self.n+2): - self.A[row][j] = -self.A[row][j] + for j in range(self.n + 2): + self.A[row][j] = -self.A[row][j] - # leave, enter in VARS defining row, col of A - # pivot tableau on the element A[row][col] which must be nonzero - # afterwards tableau normalized with positive determinant - # and updated tableau variables def pivot(self, leave, enter): n = self.n A = self.A row = self.bascobas[leave] - col = self.bascobas[enter]-n - pivelt = A[row][col] # becomes new determinant + col = self.bascobas[enter] - n + pivelt = A[row][col] negpiv = pivelt < 0 if negpiv: pivelt = -pivelt for i in range(n): if i != row: nonzero = A[i][col] != 0 - for j in range(n+2): + for j in range(n + 2): if j != col: tmp1 = A[i][j] * pivelt if nonzero: @@ -449,86 +375,99 @@ def pivot(self, leave, enter): else: tmp1 -= tmp2 A[i][j] = tmp1 // self.determinant - # row i has been dealt with, update A[i][col] safely if nonzero and not negpiv: A[i][col] = -A[i][col] - # end of for i A[row][col] = self.determinant if negpiv: self.negrow(row) - self.determinant = pivelt # by construction always positive - # update tableau variables - self.bascobas[leave] = col+n - self.whichvar[col+n] = leave + self.determinant = pivelt + self.bascobas[leave] = col + n + self.whichvar[col + n] = leave self.bascobas[enter] = row - self.whichvar[row] = enter - ###### end of pivot (leave, enter) - - def runlemke(self,*,verbose=False,lexstats=False,z0=False,silent=False): - global filehandle - # z0: printout value of z0 - # flags.maxcount = 0; - # flags.bdocupivot = 1; - # flags.binitabl = 1; - # flags.bouttabl = 0; (= verbose) - # flags.boutsol = 1; - # flags.binteract = 0; - # flags.blexstats = 0; + self.whichvar[row] = enter + def runlemke(self, *, verbose=False, lexstats=False, z0=False, silent=False): + global filehandle if silent: - filehandle = open(outfile,'w') + filehandle = open(outfile, "w") n = self.n self.pivotcount = 1 - # check if d is ok - TBC - # if (flags.binitabl) - printout ("After filltableau:") + printout("After filltableau:") printout(self) - - # z0 enters the basis to obtain lex-feasible solution enter = 0 leave, z0leave = self.lexminvar(enter) - # negate RHS - self.negcol(n+1) - # if (flags.binitabl) - if verbose: + self.negcol(n + 1) + if verbose: printout("After negcol:") printout(self) - - while True: # main loop of complementary pivoting + while True: self.testtablvars() - if z0: # printout progress of z0 - if self.bascobas[0]