diff --git a/src/lemke/bimatrix.py b/src/lemke/bimatrix.py index 268ce7b..0449b1f 100644 --- a/src/lemke/bimatrix.py +++ b/src/lemke/bimatrix.py @@ -163,7 +163,7 @@ def runLH(self, droppedlabel): if result is None: raise RuntimeError("runlemke() failed to find a solution unexpectedly.") - equilibrium = result[1: lcp.n - 1] + equilibrium = result.z[:-2] return tuple(equilibrium) def LH(self, LHstring): @@ -193,7 +193,7 @@ def runtrace(self, xprior, yprior): if result is None: raise RuntimeError("runlemke() failed to find a solution unexpectedly.") - equilibrium = result[1: lcp.n - 1] + equilibrium = result.z[:-2] return tuple(equilibrium) def trace_uniform_prior(self): diff --git a/src/lemke/lemke.py b/src/lemke/lemke.py index 5735c3d..6361263 100644 --- a/src/lemke/lemke.py +++ b/src/lemke/lemke.py @@ -3,6 +3,7 @@ import fractions import math # gcd import sys +from dataclasses import dataclass import click @@ -99,7 +100,7 @@ def __init__(self, Mqd): 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 @@ -176,21 +177,6 @@ def vartoa(self, v): # variable as as string w1..wn or z0..zn else: return "z" + str(v) - def createsol(self): # get solution from current tableau - n = self.n - 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 assertbasic(self, v, info): # assert that v is basic if self.bascobas[v] >= self.n: raise RuntimeError( @@ -355,35 +341,13 @@ def pivot(self, leave, enter): class RayTermination(Exception): def __init__(self, enter, tableau): - tableau.createsol() self.tableau = tableau + self.enter = tableau.vartoa(enter) super().__init__( - "Ray termination when trying to enter " + tableau.vartoa(enter) + "Ray termination when trying to enter " + self.enter ) -def outsol(tableau): # string giving solution, after createsol() - # printout in columns to check complementarity - n = tableau.n - sol = columnprint.columnprint(n + 2) - sol.sprint("basis=") - for i in range(n + 1): - if tableau.bascobas[i] < n: # Z(i) is a basic variable - s = tableau.vartoa(i) - elif i > 0 and tableau.bascobas[n + i] < n: # W(i) is a basic variable - s = tableau.vartoa(n + i) - else: - s = " " - sol.sprint(s) - sol.sprint("z=") - for i in range(2 * n + 1): - sol.sprint(str(tableau.solution[i])) - if i == n: # new line since printouting slack vars w next - sol.sprint("w=") - sol.sprint("") # no W(0) - return str(sol) - - # output statistics of minimum ratio test def outstatistics(tableau): n = tableau.n @@ -415,8 +379,8 @@ def on_start(self, lcp, tableau): pass def on_negcol(self, tableau): pass def on_pivot_start(self, tableau, leave, enter): pass def on_pivot_end(self, tableau): pass - def on_done(self, tableau): pass - def on_ray_termination(self, tableau, message): pass + def on_done(self, tableau, result): pass + def on_ray_termination(self, tableau, result, message): pass class PrintingCallback(LemkeCallback): @@ -473,7 +437,7 @@ def on_pivot_end(self, tableau): if self.verbose: self.printout(tableau) - def on_done(self, tableau): + def on_done(self, tableau, result): if self.z0: self.printout(f"pivot count = {tableau.pivotcount + 1}, z0 = 0.0") @@ -482,22 +446,118 @@ def on_done(self, tableau): self.printout(tableau) # if (flags.boutsol) - self.printout(outsol(tableau)) + self.printout(result) if self.lexstats: # output statistics of minimum ratio test self.printout(outstatistics(tableau)) - def on_ray_termination(self, tableau, message): + def on_ray_termination(self, tableau, result, message): self.printout(message) self.printout(tableau) - self.printout("Current basis not an LCP solution:") - self.printout(outsol(tableau)) + self.printout(result) + + +@dataclass(frozen=True) +class LcpResult: + success: bool + num_pivots: int + basis: frozenset[str] # e.g. {'w1', 'z2', ...} + z0: fractions.Fraction + z: tuple[fractions.Fraction, ...] # (z1, ..., zn) + w: tuple[fractions.Fraction, ...] # (w1, ..., wn) + ray_entering_variable: str | None + + def __str__(self): + pivot_word = "pivot" if self.num_pivots == 1 else "pivots" + + if self.success: + status = ( + f"Process finished successfully after {self.num_pivots} {pivot_word}.\n" + "Solution found:\n" + ) + else: + status = ( + f"Terminated on a secondary ray after {self.num_pivots} {pivot_word}, " + f"when trying to enter {self.ray_entering_variable}.\n" + "Current basis not an LCP solution:\n" + ) + + # printout in columns to check complementarity + n = len(self.w) + + sol = columnprint.columnprint(n + 2) + + sol.sprint("basis=") + # align basis elements with corresponding columns + basis_by_row = {int(b[1:]): b for b in self.basis} + for i in range(n + 1): + sol.sprint(basis_by_row.get(i, " ")) + + sol.sprint("z=") + sol.sprint(str(self.z0)) + for el in self.z: + sol.sprint(str(el)) + + sol.sprint("w=") + sol.sprint("") # no W(0) + for el in self.w: + sol.sprint(str(el)) + + return status + str(sol) + + +def result_from_tableau( + tableau: tableau, + success: bool, + ray_entering_variable: str | None = None, +) -> LcpResult: + n = tableau.n + basis = set() + + # [z0, z1, ..., zn, w1, ..., wn] + solution = [fractions.Fraction(0) for _ in range(2 * n + 1)] + + for i in range(2 * n + 1): + row = tableau.bascobas[i] + if row < n: # i is a basic variable + num = tableau.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 *= tableau.scalefactor[i] + solution[i] = fractions.Fraction( + num, + tableau.determinant * tableau.scalefactor[n + 1] + ) + basis.add(tableau.vartoa(i)) + + return LcpResult( + success=success, + num_pivots=tableau.pivotcount, + basis=frozenset(basis), + z0=solution[0], + z=tuple(solution[1:n + 1]), + w=tuple(solution[n + 1:]), + ray_entering_variable=ray_entering_variable, + ) def runlemke(*, lcp, callback=None): callback = callback or LemkeCallback() + # trivial case (q >= 0) + if all(element >= 0 for element in lcp.q): + return LcpResult( + success=True, + num_pivots=0, + basis=frozenset(f"w{i + 1}" for i in range(lcp.n)), + z0=fractions.Fraction(0), + z=(fractions.Fraction(0),) * lcp.n, + w=tuple(lcp.q), + ray_entering_variable=None, + ) + try: tabl = tableau(lcp) @@ -534,13 +594,18 @@ def runlemke(*, lcp, callback=None): leave, z0leave = tabl.lexminvar(enter) tabl.pivotcount += 1 - tabl.createsol() - callback.on_done(tableau=tabl) + result = result_from_tableau(tabl, True) + callback.on_done(tableau=tabl, result=result) - return tabl.solution + return result except RayTermination as e: - callback.on_ray_termination(message=str(e), tableau=e.tableau) - return None + result = result_from_tableau( + tableau=e.tableau, + success=False, + ray_entering_variable=e.enter, + ) + callback.on_ray_termination(message=str(e), result=result, tableau=e.tableau) + return result @click.command( @@ -575,7 +640,7 @@ def main(verbose, z0, lcpfilename): callback=PrintingCallback(stream=sys.stdout, verbose=verbose, z0=z0), ) - if result is None: + if not result.success: sys.exit(1) diff --git a/tests/sequence_form_helper.py b/tests/sequence_form_helper.py index 41c683f..b1845a2 100644 --- a/tests/sequence_form_helper.py +++ b/tests/sequence_form_helper.py @@ -223,10 +223,10 @@ def solve_via_sequence_form(game): d = [Fr(1) for _ in range(len(q))] lcp_instance = lcp_from_data(M, q, d) - sol = runlemke(lcp=lcp_instance) + result = runlemke(lcp=lcp_instance) # realization plans - x_y = sol[1:(ns1 + ns2 + 1)] + x_y = result.z x = x_y[:ns1] y = x_y[ns1:] diff --git a/tests/test_lcp.py b/tests/test_lcp.py index 8b2d4e4..7effcad 100644 --- a/tests/test_lcp.py +++ b/tests/test_lcp.py @@ -26,7 +26,9 @@ class LCPTestCase: """Defines data for one LCP test case for Lemke's algorithm.""" factory: Callable[[], lcp] - expected: list[Fr] | None = None + expected_z: list[Fr] | None = None + expected_w: list[Fr] | None = None + expected_basis: set[str] | None = None tol: Fr = Fr(0) @@ -37,7 +39,9 @@ class LCPTestCase: pytest.param( LCPTestCase( factory=lambda: lcp.from_file(FIXTURES_DIR / "trivial_q_pos_M_arbitrary"), - expected=[Fr(0), Fr(0), Fr(0), Fr(3), Fr(1)], + expected_z=[Fr(0), Fr(0)], + expected_w=[Fr(3), Fr(1)], + expected_basis={"w1", "w2"}, ), id="trivial_q_pos_M_arbitrary", ), @@ -45,7 +49,9 @@ class LCPTestCase: pytest.param( LCPTestCase( factory=lambda: lcp.from_file(FIXTURES_DIR / "trivial_q_zero_M_identity"), - expected=[Fr(0), Fr(0), Fr(0), Fr(0), Fr(0)], + expected_z=[Fr(0), Fr(0)], + expected_w=[Fr(0), Fr(0)], + expected_basis={"w1", "w2"}, ), id="trivial_q_zero_M_identity", ), @@ -53,7 +59,9 @@ class LCPTestCase: pytest.param( LCPTestCase( factory=lambda: lcp.from_file(FIXTURES_DIR / "trivial_q_zero_M_zero"), - expected=[Fr(0), Fr(0), Fr(0), Fr(0), Fr(0)], + expected_z=[Fr(0), Fr(0)], + expected_w=[Fr(0), Fr(0)], + expected_basis={"w1", "w2"}, ), id="trivial_q_zero_M_zero", ), @@ -67,7 +75,9 @@ class LCPTestCase: pytest.param( LCPTestCase( factory=lambda: lcp.from_file(FIXTURES_DIR / "non_degenerate_book_ex_3x3"), - expected=[Fr(0), Fr(0), Fr(1), Fr(3), Fr(2), Fr(0), Fr(0)], + expected_z=[Fr(0), Fr(1), Fr(3)], + expected_w=[Fr(2), Fr(0), Fr(0)], + expected_basis={"w1", "z2", "z3"}, ), id="non_degenerate_book_ex_3x3", ), @@ -76,7 +86,9 @@ class LCPTestCase: pytest.param( LCPTestCase( factory=lambda: lcp.from_file(FIXTURES_DIR / "non_degenerate_book_ex_4x4"), - expected=[Fr(0), Fr(0), Fr(1, 2), Fr(0), Fr(0), Fr(1, 2), Fr(0), Fr(11, 2), Fr(4)], + expected_z=[Fr(0), Fr(1, 2), Fr(0), Fr(0)], + expected_w=[Fr(1, 2), Fr(0), Fr(11, 2), Fr(4)], + expected_basis={"w1", "z2", "w3", "w4"}, ), id="non_degenerate_book_ex_4x4", ), @@ -86,7 +98,9 @@ class LCPTestCase: pytest.param( LCPTestCase( factory=lambda: lcp.from_file(FIXTURES_DIR / "non_degenerate_1x1"), - expected=[Fr(0), Fr(329, 20), Fr(0)], + expected_z=[Fr(329, 20)], + expected_w=[Fr(0)], + expected_basis={"z1"}, ), id="non_degenerate_1x1", ), @@ -94,7 +108,9 @@ class LCPTestCase: pytest.param( LCPTestCase( factory=lambda: lcp.from_file(FIXTURES_DIR / "non_degenerate_2x2"), - expected=[Fr(0), Fr(2), Fr(1), Fr(0), Fr(0)], + expected_z=[Fr(2), Fr(1)], + expected_w=[Fr(0), Fr(0)], + expected_basis={"z1", "z2"}, ), id="non_degenerate_2x2", ), @@ -104,7 +120,9 @@ class LCPTestCase: pytest.param( LCPTestCase( factory=lambda: lcp.from_file(FIXTURES_DIR / "non_degenerate_M_identity"), - expected=[Fr(0), Fr(51, 10), Fr(0), Fr(0), Fr(8), Fr(0), Fr(2, 7), Fr(10), Fr(0)], + expected_z=[Fr(51, 10), Fr(0), Fr(0), Fr(8)], + expected_w=[Fr(0), Fr(2, 7), Fr(10), Fr(0)], + expected_basis={"z1", "w2", "w3", "z4"}, ), id="non_degenerate_M_identity", ), @@ -119,7 +137,9 @@ class LCPTestCase: pytest.param( LCPTestCase( factory=lambda: lcp.from_file(FIXTURES_DIR / "degenerate_tie_in_initial_lexmin"), - expected=[Fr(0), Fr(0), Fr(1), Fr(0), Fr(0)], + expected_z=[Fr(0), Fr(1)], + expected_w=[Fr(0), Fr(0)], + expected_basis={"w1", "z2"}, ), id="degenerate_tie_in_initial_lexmin", ), @@ -127,7 +147,9 @@ class LCPTestCase: pytest.param( LCPTestCase( factory=lambda: lcp.from_file(FIXTURES_DIR / "degenerate_tie_in_noninitial_lexmin"), - expected=[Fr(0), Fr(0), Fr(1), Fr(0), Fr(0)], + expected_z=[Fr(0), Fr(1)], + expected_w=[Fr(0), Fr(0)], + expected_basis={"w1", "z2"}, ), id="degenerate_tie_in_noninitial_lexmin", ), @@ -135,7 +157,9 @@ class LCPTestCase: pytest.param( LCPTestCase( factory=lambda: lcp.from_file(FIXTURES_DIR / "degenerate_tie_in_several_lexmins"), - expected=[Fr(0), Fr(0), Fr(1), Fr(0), Fr(1), Fr(0), Fr(0)], + expected_z=[Fr(0), Fr(1), Fr(0)], + expected_w=[Fr(1), Fr(0), Fr(0)], + expected_basis={"w1", "z2", "w3"}, ), id="degenerate_tie_in_several_lexmins", ), @@ -143,7 +167,7 @@ class LCPTestCase: SUCCESS_CASES = [] -# SUCCESS_CASES += TRIVIAL_CASES +SUCCESS_CASES += TRIVIAL_CASES SUCCESS_CASES += NON_DEGENERATE_CASES SUCCESS_CASES += DEGENERATE_CASES @@ -156,17 +180,25 @@ def test_with_expected_results(test_case: LCPTestCase, subtests): """ lcp_instance = test_case.factory() - sol = runlemke(lcp=lcp_instance) - n = lcp_instance.n + result = runlemke(lcp=lcp_instance) + + with subtests.test("Solution status"): + assert result.success - with subtests.test("Solution length"): - assert len(sol) == 2 * n + 1 + with subtests.test("z0 value"): + assert result.z0 == Fr(0) - for i, val in enumerate(sol): - label = f"z{i}" if i <= n else f"w{i - n}" - expected_val = test_case.expected[i] + with subtests.test("Basis value"): + assert result.basis == test_case.expected_basis + + for i, val in enumerate(result.z): + expected_val = test_case.expected_z[i] + with subtests.test(f"z{i + 1} value"): + assert abs(val - expected_val) <= test_case.tol - with subtests.test(f"{label} value"): + for i, val in enumerate(result.w): + expected_val = test_case.expected_w[i] + with subtests.test(f"w{i + 1} value"): assert abs(val - expected_val) <= test_case.tol @@ -182,17 +214,35 @@ def test_with_lcp_conditions(test_case: LCPTestCase, subtests): """ lcp_instance = test_case.factory() - sol = runlemke(lcp=lcp_instance) + result = runlemke(lcp=lcp_instance) - # solution format: [z0, z1..zn, w1..wn] n = lcp_instance.n - z0 = sol[0] - z = sol[1:n + 1] - w = sol[n + 1:] + z0 = result.z0 + z = result.z + w = result.w + + with subtests.test("Solution status"): + assert result.success with subtests.test("z0 = 0"): assert z0 == Fr(0) + for i in range(n): + with subtests.test(index=i): + occurrences = sum( + 1 for var in result.basis + if int(var[1:]) == i + ) + + if i == 0: + assert occurrences <= 1, ( + f"Index {i} appears {occurrences} times in basis" + ) + else: + assert occurrences == 1, ( + f"Index {i} appears {occurrences} times in basis" + ) + for i, val in enumerate(z): with subtests.test(f"z{i + 1} nonnegativity"): assert val >= 0 @@ -239,4 +289,6 @@ def test_failure(test_case: LCPTestCase): """ lcp_instance = test_case.factory() - assert runlemke(lcp=lcp_instance) is None + result = runlemke(lcp=lcp_instance) + + assert result.success is False