From b0e8b8689be6a45131078eee1c06a994225b09fe Mon Sep 17 00:00:00 2001
From: nataliemes <nataliemsbl@gmail.com>
Date: Tue, 9 Jun 2026 22:22:32 +0400
Subject: [PATCH 1/8] Separate runlemke() from tableau class

---
 src/lemke/bimatrix.py         |  15 ++--
 src/lemke/lemke.py            | 127 ++++++++++++++++++----------------
 tests/sequence_form_helper.py |   5 +-
 tests/test_lcp.py             |  16 ++---
 4 files changed, 76 insertions(+), 87 deletions(-)

diff --git a/src/lemke/bimatrix.py b/src/lemke/bimatrix.py
index 241e34c..5c34d49 100644
--- a/src/lemke/bimatrix.py
+++ b/src/lemke/bimatrix.py
@@ -251,10 +251,9 @@ def createLCP(self):
     def runLH(self, droppedlabel):
         lcp = self.createLCP()
         lcp.d[droppedlabel - 1] = 0  # subsidize this label
-        tabl = lemke.tableau(lcp)
         # tabl.runlemke(verbose=True, lexstats=True, z0=gz0)
-        tabl.runlemke(silent=True)
-        return tuple(getequil(tabl))
+        equilibrium = lemke.runlemke(lcp=lcp, silent=True)[1: lcp.n - 1]
+        return tuple(equilibrium)
 
     def LH(self, LHstring):
         if LHstring == "":
@@ -280,9 +279,8 @@ def runtrace(self, xprior, yprior):
         Ay = self.A.negmatrix @ yprior
         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))
+        equilibrium = lemke.runlemke(lcp=lcp, silent=True)[1: lcp.n - 1]
+        return tuple(equilibrium)
 
     def tracing(self, trace):
         if trace < 0:
@@ -337,11 +335,6 @@ def uniform(n):
     return np.array([fractions.Fraction(1, n) for _ 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]]) + ")"
diff --git a/src/lemke/lemke.py b/src/lemke/lemke.py
index 4aad705..acf2193 100644
--- a/src/lemke/lemke.py
+++ b/src/lemke/lemke.py
@@ -468,68 +468,74 @@ def pivot(self, leave, 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;
-
-        if silent:
-            filehandle = open(outfile, "w")  # noqa: SIM115
-        n = self.n
-        self.pivotcount = 1
-        # check if d is ok - TBC
-        # if (flags.binitabl)
-        printout("After filltableau:")
-        printout(self)
+    #  end of class tableau
+
 
-        # 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)
+def runlemke(*, lcp, 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;
+
+    tabl = tableau(lcp)
+
+    if silent:
+        filehandle = open(outfile, "w")  # noqa: SIM115
+    n = tabl.n
+    tabl.pivotcount = 1
+    # check if d is ok - TBC
+    # if (flags.binitabl)
+    printout("After filltableau:")
+    printout(tabl)
+
+    # z0 enters the basis to obtain lex-feasible solution
+    enter = 0
+    leave, z0leave = tabl.lexminvar(enter)
+    # negate RHS
+    tabl.negcol(n + 1)
+    # if (flags.binitabl)
+    if verbose:
+        printout("After negcol:")
+        printout(tabl)
+
+    while True:  # main loop of complementary pivoting
+        tabl.testtablvars()
+        if z0:  # printout progress of z0
+            if tabl.bascobas[0] < n:  # z0 is basic
+                printout(
+                    "step,z0=",
+                    tabl.pivotcount, tabl.A[tabl.bascobas[0]][n + 1] / tabl.determinant
+                )
+            else:
+                printout("step,z0=", tabl.pivotcount, 0.0)
+        # if (flags.bdocupivot)
+        tabl.docupivot(leave, enter)
+        tabl.pivot(leave, enter)
+        if z0leave:
+            if z0:
+                printout("step,z0=", tabl.pivotcount + 1, 0.0)
+            break
         if verbose:
-            printout("After negcol:")
-            printout(self)
+            printout(tabl)
+        enter = tabl.complement(leave)
+        leave, z0leave = tabl.lexminvar(enter)
+        tabl.pivotcount += 1
 
-        while True:  # main loop of complementary pivoting
-            self.testtablvars()
-            if z0:  # printout progress of z0
-                if self.bascobas[0] < n:  # z0 is basic
-                    printout(
-                        "step,z0=",
-                        self.pivotcount, self.A[self.bascobas[0]][n + 1] / self.determinant
-                    )
-                else:
-                    printout("step,z0=", self.pivotcount, 0.0)
-            # if (flags.bdocupivot)
-            self.docupivot(leave, enter)
-            self.pivot(leave, enter)
-            if z0leave:
-                if z0:
-                    printout("step,z0=", self.pivotcount + 1, 0.0)
-                break
-            if verbose:
-                printout(self)
-            enter = self.complement(leave)
-            leave, z0leave = self.lexminvar(enter)
-            self.pivotcount += 1
-
-        # if (flags.binitabl)
-        printout("Final tableau:")
-        printout(self)
-        # if (flags.boutsol)
-        self.createsol()
-        printout(self.outsol())
-        if lexstats:
-            self.outstatistics()
-    #  end of class tableau
+    # if (flags.binitabl)
+    printout("Final tableau:")
+    printout(tabl)
+    # if (flags.boutsol)
+    tabl.createsol()
+    printout(tabl.outsol())
+    if lexstats:
+        tabl.outstatistics()
+
+    return tabl.solution
 
 
 def main():
@@ -539,8 +545,7 @@ def main():
     m = lcp(lcpfilename)
     printout(m)
     printout("==================================")
-    tabl = tableau(m)
-    tabl.runlemke(verbose=verbose, z0=z0, silent=silent)
+    runlemke(lcp=m, verbose=verbose, z0=z0, silent=silent)
 
 
 if __name__ == "__main__":
diff --git a/tests/sequence_form_helper.py b/tests/sequence_form_helper.py
index a9a0cd5..d9cd8ea 100644
--- a/tests/sequence_form_helper.py
+++ b/tests/sequence_form_helper.py
@@ -1,8 +1,7 @@
 from collections import defaultdict
 from fractions import Fraction as Fr
 
-from lemke.lemke import lcp
-from tests.test_lcp import lemke_solver
+from lemke.lemke import lcp, runlemke
 
 EMPTY: tuple[()] = ()  # empty sequence
 
@@ -221,7 +220,7 @@ def solve_via_sequence_form(game):
     d = [Fr(1) for _ in range(len(q))]
 
     lcp_instance = lcp_from_data(M, q, d)
-    sol = lemke_solver(lcp_instance)
+    sol = runlemke(lcp=lcp_instance)
 
     # realization plans
     x_y = sol[1:(ns1 + ns2 + 1)]
diff --git a/tests/test_lcp.py b/tests/test_lcp.py
index 4198b85..d856f93 100644
--- a/tests/test_lcp.py
+++ b/tests/test_lcp.py
@@ -17,19 +17,11 @@
 
 import pytest
 
-from lemke.lemke import lcp, tableau
+from lemke.lemke import lcp, runlemke
 
 FIXTURES_DIR = Path("tests/fixtures/lcp")
 
 
-def lemke_solver(lcp_instance: lcp) -> list[Fr]:
-    """Runs Lemke's algorithm on the given LCP and returns the solution."""
-
-    tabl = tableau(lcp_instance)
-    tabl.runlemke(verbose=False, z0=False, silent=False)
-    return tabl.solution
-
-
 @dataclass
 class LCPTestCase:
     """Defines data for one LCP test case for Lemke's algorithm."""
@@ -164,7 +156,7 @@ def test_with_expected_results(test_case: LCPTestCase, subtests):
     """
 
     lcp_instance = test_case.factory()
-    sol = lemke_solver(lcp_instance)
+    sol = runlemke(lcp=lcp_instance)
     n = lcp_instance.n
 
     with subtests.test("Solution length"):
@@ -190,7 +182,7 @@ def test_with_lcp_conditions(test_case: LCPTestCase, subtests):
     """
 
     lcp_instance = test_case.factory()
-    sol = lemke_solver(lcp_instance)
+    sol = runlemke(lcp=lcp_instance)
 
     # solution format: [z0, z1..zn, w1..wn]
     n = lcp_instance.n
@@ -248,6 +240,6 @@ def test_failure(test_case: LCPTestCase):
     lcp_instance = test_case.factory()
 
     with pytest.raises(SystemExit) as exc_info:
-        lemke_solver(lcp_instance)
+        runlemke(lcp=lcp_instance)
 
     assert exc_info.value.code == 1

From 3de14b7c65e66ee1563d76e94c8a6ff58809a67b Mon Sep 17 00:00:00 2001
From: nataliemes <nataliemsbl@gmail.com>
Date: Tue, 16 Jun 2026 13:44:25 +0400
Subject: [PATCH 2/8] Remove printout() calls outside runlemke()

---
 src/lemke/lemke.py | 197 +++++++++++++++++++++++----------------------
 tests/test_lcp.py  |   5 +-
 2 files changed, 103 insertions(+), 99 deletions(-)

diff --git a/src/lemke/lemke.py b/src/lemke/lemke.py
index acf2193..1a3c557 100644
--- a/src/lemke/lemke.py
+++ b/src/lemke/lemke.py
@@ -40,7 +40,7 @@ def processArguments():
             lcpfilename = s
             outfile = s + ".out"
     if showhelp:
-        printout(helpstring)
+        print(helpstring)
         exit(0)
     return
 
@@ -70,10 +70,10 @@ def __init__(self, arg):
             # 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]))
-                exit(1)
+                raise ValueError(
+                    f"lcp file {filename!r} must start with 'n=' lcpdim, e.g. 'n= 5', "
+                    f"not {words[0]!r}"
+                )
             n = int(words[1])
             self.n = n
             # self.M = np.zeros( (n,n), dtype=fractions.Fraction)
@@ -87,10 +87,11 @@ def __init__(self, arg):
             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)
-                exit(1)
+                raise ValueError(
+                    f"in lcp file {filename!r}: "
+                    f"n={n}, need keywords 'M=' 'q=' 'd=' and "
+                    f"n*n + n + n = {needfracs} fractions, got {len(words) - 5}"
+                )
             k = 2  # index in words
             while k < len(words):
                 if words[k] == "M=":
@@ -106,9 +107,10 @@ def __init__(self, arg):
                     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)
+                    raise ValueError(
+                        f"in lcp file {filename!r}: expected one of 'M=' 'q=' 'd=', "
+                        f"got {words[k]!r}"
+                    )
             return
 
     def __str__(self):
@@ -266,54 +268,42 @@ def outsol(self):  # string giving solution, after createsol()
 
     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
+            raise RuntimeError(
+                f"({info}) Cobasic variable {self.vartoa(v)} should be basic"
+            )
 
     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
+            raise RuntimeError(
+                f"({info}) Basic variable {self.vartoa(v)} should be cobasic"
+            )
 
     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
-
-    def raytermination(self, enter):
-        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)
+        return s
 
-    def testtablvars(self):  # msg only if error, continue
+    def testtablvars(self):
         n = self.n
         for i in range(2 * n + 1):
             if self.bascobas[self.whichvar[i]] != i:
+                message = ""
                 # 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(
+                        message += f"First problem for j={j}:\n"
+                    message += (
                         f"j={j} self.bascobas[j]={self.bascobas[j]} "
-                        f"self.whichvar[j]={self.whichvar[j]}"
+                        f"self.whichvar[j]={self.whichvar[j]}\n"
                     )
-                break
+                raise RuntimeError(f"testtablvars() failed:\n{message}")
 
     def complement(self, v):  # Z(i),W(i) are complements
         n = self.n
         if v == 0:
-            printout("Attempt to find complement of z0")
-            exit(1)
+            raise RuntimeError("Attempt to find complement of z0")
         if v > n:
             return v - n
         else:
@@ -342,7 +332,7 @@ def outstatistics(self):
                     stats.sprint(str(x / 10.0))
                 else:
                     stats.sprint("-")
-        printout(stats)
+        return stats
 
     # returns leave,z0leave
     # leave = leaving variable in VARS, given by lexmin row,
@@ -361,7 +351,7 @@ def lexminvar(self, enter):
             if A[i][col] > 0:
                 leavecand.append(i)
         if not leavecand:
-            self.raytermination(enter)
+            raise RayTermination(enter, self)
         if len(leavecand) == 1:  # single positive entering value
             z0leave = self.bascobas[0] == leavecand[0]
         # omitted from statistics: only one possible row
@@ -377,8 +367,7 @@ def lexminvar(self, enter):
         j = 0  # going through j = 0..n
         while len(leavecand) > 1:
             if j > n:  # impossible, perturbed RHS should have full rank
-                printout("lex-minratio test failed")
-                exit(1)
+                raise RuntimeError("lex-minratio test failed")
             self.lextested[j] += 1
             self.lexcomparisons[j] += len(leavecand)
             testcol = n + 1 if j == 0 else self.bascobas[n + j] - n
@@ -471,6 +460,15 @@ def pivot(self, leave, enter):
     #  end of class tableau
 
 
+class RayTermination(Exception):
+    def __init__(self, enter, tableau):
+        tableau.createsol()
+        self.tableau = tableau
+        super().__init__(
+            "Ray termination when trying to enter " + tableau.vartoa(enter)
+        )
+
+
 def runlemke(*, lcp, verbose=False, lexstats=False, z0=False, silent=False):
     global filehandle
     # z0: printout value of z0
@@ -482,69 +480,78 @@ def runlemke(*, lcp, verbose=False, lexstats=False, z0=False, silent=False):
     # flags.binteract  = 0;
     # flags.blexstats  = 0;
 
-    tabl = tableau(lcp)
-
-    if silent:
-        filehandle = open(outfile, "w")  # noqa: SIM115
-    n = tabl.n
-    tabl.pivotcount = 1
-    # check if d is ok - TBC
-    # if (flags.binitabl)
-    printout("After filltableau:")
-    printout(tabl)
-
-    # z0 enters the basis to obtain lex-feasible solution
-    enter = 0
-    leave, z0leave = tabl.lexminvar(enter)
-    # negate RHS
-    tabl.negcol(n + 1)
-    # if (flags.binitabl)
-    if verbose:
-        printout("After negcol:")
+    try:
+        tabl = tableau(lcp)
+
+        printout(f"verbose={verbose} lcpfilename={lcpfilename} silent={silent} z0={z0}")
+        printout(lcp)
+        printout("==================================")
+
+        if silent:
+            filehandle = open(outfile, "w")  # noqa: SIM115
+        n = tabl.n
+        tabl.pivotcount = 1
+        # check if d is ok - TBC
+        # if (flags.binitabl)
+        printout("After filltableau:")
         printout(tabl)
 
-    while True:  # main loop of complementary pivoting
-        tabl.testtablvars()
-        if z0:  # printout progress of z0
-            if tabl.bascobas[0] < n:  # z0 is basic
-                printout(
-                    "step,z0=",
-                    tabl.pivotcount, tabl.A[tabl.bascobas[0]][n + 1] / tabl.determinant
-                )
-            else:
-                printout("step,z0=", tabl.pivotcount, 0.0)
-        # if (flags.bdocupivot)
-        tabl.docupivot(leave, enter)
-        tabl.pivot(leave, enter)
-        if z0leave:
-            if z0:
-                printout("step,z0=", tabl.pivotcount + 1, 0.0)
-            break
+        # z0 enters the basis to obtain lex-feasible solution
+        enter = 0
+        leave, z0leave = tabl.lexminvar(enter)
+        # negate RHS
+        tabl.negcol(n + 1)
+        # if (flags.binitabl)
         if verbose:
+            printout("After negcol:")
             printout(tabl)
-        enter = tabl.complement(leave)
-        leave, z0leave = tabl.lexminvar(enter)
-        tabl.pivotcount += 1
 
-    # if (flags.binitabl)
-    printout("Final tableau:")
-    printout(tabl)
-    # if (flags.boutsol)
-    tabl.createsol()
-    printout(tabl.outsol())
-    if lexstats:
-        tabl.outstatistics()
-
-    return tabl.solution
+        while True:  # main loop of complementary pivoting
+            tabl.testtablvars()
+            if z0:  # printout progress of z0
+                if tabl.bascobas[0] < n:  # z0 is basic
+                    printout(
+                        "step,z0=",
+                        tabl.pivotcount, tabl.A[tabl.bascobas[0]][n + 1] / tabl.determinant
+                    )
+                else:
+                    printout("step,z0=", tabl.pivotcount, 0.0)
+            # if (flags.bdocupivot)
+            printout(tabl.docupivot(leave, enter))
+            tabl.pivot(leave, enter)
+            if z0leave:
+                if z0:
+                    printout("step,z0=", tabl.pivotcount + 1, 0.0)
+                break
+            if verbose:
+                printout(tabl)
+            enter = tabl.complement(leave)
+            leave, z0leave = tabl.lexminvar(enter)
+            tabl.pivotcount += 1
+
+        # if (flags.binitabl)
+        printout("Final tableau:")
+        printout(tabl)
+        # if (flags.boutsol)
+        tabl.createsol()
+        printout(tabl.outsol())
+        if lexstats:
+            printout(tabl.outstatistics())
+
+        return tabl.solution
+    except RayTermination as e:
+        printout(str(e))
+        printout(e.tableau)
+        printout("Current basis not an LCP solution:")
+        printout(e.tableau.outsol())
+        return None
 
 
 def main():
     processArguments()
-    printout(f"verbose={verbose} lcpfilename={lcpfilename} silent={silent} z0={z0}")
-    # printout (f"{verbose}= {lcpfilename}= {silent}= {z0}=")
+
     m = lcp(lcpfilename)
-    printout(m)
-    printout("==================================")
+
     runlemke(lcp=m, verbose=verbose, z0=z0, silent=silent)
 
 
diff --git a/tests/test_lcp.py b/tests/test_lcp.py
index d856f93..4f225b9 100644
--- a/tests/test_lcp.py
+++ b/tests/test_lcp.py
@@ -239,7 +239,4 @@ def test_failure(test_case: LCPTestCase):
     """
     lcp_instance = test_case.factory()
 
-    with pytest.raises(SystemExit) as exc_info:
-        runlemke(lcp=lcp_instance)
-
-    assert exc_info.value.code == 1
+    assert runlemke(lcp=lcp_instance) is None

From 31676360f7ee315ca2ec9f6783bf4fcca713310c Mon Sep 17 00:00:00 2001
From: nataliemes <nataliemsbl@gmail.com>
Date: Tue, 16 Jun 2026 13:53:24 +0400
Subject: [PATCH 3/8] Separate docupivot(), outsol(), outstatistics() from
 tableau class

---
 src/lemke/lemke.py | 115 +++++++++++++++++++++++----------------------
 1 file changed, 58 insertions(+), 57 deletions(-)

diff --git a/src/lemke/lemke.py b/src/lemke/lemke.py
index 1a3c557..cee3427 100644
--- a/src/lemke/lemke.py
+++ b/src/lemke/lemke.py
@@ -245,27 +245,6 @@ def createsol(self):  # get solution from current tableau
             else:  # i is nonbasic
                 self.solution[i] = fractions.Fraction(0)
 
-    def outsol(self):  # string giving solution, after createsol()
-        # printout in columns to check complementarity
-        n = self.n
-        sol = columnprint.columnprint(n + 2)
-        sol.sprint("basis=")
-        for i in range(n + 1):
-            if self.bascobas[i] < n:  # Z(i) is a basic variable
-                s = self.vartoa(i)
-            elif 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):
-            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:
             raise RuntimeError(
@@ -278,13 +257,6 @@ def assertcobasic(self, v, info):  # assert that v is cobasic
                 f"({info}) Basic variable {self.vartoa(v)} should be cobasic"
             )
 
-    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)
-        return s
-
     def testtablvars(self):
         n = self.n
         for i in range(2 * n + 1):
@@ -309,31 +281,6 @@ def complement(self, v):  # Z(i),W(i) are complements
         else:
             return v + n
 
-    # output statistics of minimum ratio test
-    def outstatistics(self):
-        n = self.n
-        lext = self.lextested
-        stats = columnprint.columnprint(n + 2)
-        stats.makeLeft(0)
-        stats.sprint("lex-column")
-        for i in range(n + 1):
-            stats.iprint(i)
-        stats.sprint("times tested")
-        for i in range(n + 1):
-            stats.iprint(lext[i])
-        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))
-            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))
-                else:
-                    stats.sprint("-")
-        return stats
-
     # returns leave,z0leave
     # leave = leaving variable in VARS, given by lexmin row,
     # when enter in VARS is entering variable
@@ -469,6 +416,58 @@ def __init__(self, enter, tableau):
         )
 
 
+def docupivot(leave, enter):  # leave, enter in VARS
+    return f"leaving: {leave.ljust(5)} entering: {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
+    lext = tableau.lextested
+    stats = columnprint.columnprint(n + 2)
+    stats.makeLeft(0)
+    stats.sprint("lex-column")
+    for i in range(n + 1):
+        stats.iprint(i)
+    stats.sprint("times tested")
+    for i in range(n + 1):
+        stats.iprint(lext[i])
+    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 / tableau.pivotcount))
+        stats.sprint("avg comparisons")
+        for i in range(n + 1):
+            if lext[i] > 0:
+                x = round(tableau.lexcomparisons[i] * 10 / lext[0])
+                stats.sprint(str(x / 10.0))
+            else:
+                stats.sprint("-")
+    return stats
+
+
 def runlemke(*, lcp, verbose=False, lexstats=False, z0=False, silent=False):
     global filehandle
     # z0: printout value of z0
@@ -517,7 +516,9 @@ def runlemke(*, lcp, verbose=False, lexstats=False, z0=False, silent=False):
                 else:
                     printout("step,z0=", tabl.pivotcount, 0.0)
             # if (flags.bdocupivot)
-            printout(tabl.docupivot(leave, enter))
+            tabl.assertbasic(leave, "docupivot")
+            tabl.assertcobasic(enter, "docupivot")
+            printout(docupivot(tabl.vartoa(leave), tabl.vartoa(enter)))
             tabl.pivot(leave, enter)
             if z0leave:
                 if z0:
@@ -534,16 +535,16 @@ def runlemke(*, lcp, verbose=False, lexstats=False, z0=False, silent=False):
         printout(tabl)
         # if (flags.boutsol)
         tabl.createsol()
-        printout(tabl.outsol())
+        printout(outsol(tabl))
         if lexstats:
-            printout(tabl.outstatistics())
+            printout(outstatistics(tabl))
 
         return tabl.solution
     except RayTermination as e:
         printout(str(e))
         printout(e.tableau)
         printout("Current basis not an LCP solution:")
-        printout(e.tableau.outsol())
+        printout(outsol(e.tableau))
         return None
 
 

From 945bc47850ceaa2c16d11b645d980d0badec5b69 Mon Sep 17 00:00:00 2001
From: nataliemes <nataliemsbl@gmail.com>
Date: Tue, 16 Jun 2026 14:40:46 +0400
Subject: [PATCH 4/8] Add callback functions

---
 src/lemke/bimatrix.py |   4 +-
 src/lemke/lemke.py    | 152 ++++++++++++++++++++++++++++--------------
 2 files changed, 105 insertions(+), 51 deletions(-)

diff --git a/src/lemke/bimatrix.py b/src/lemke/bimatrix.py
index 5c34d49..e7cdb20 100644
--- a/src/lemke/bimatrix.py
+++ b/src/lemke/bimatrix.py
@@ -252,7 +252,7 @@ def runLH(self, droppedlabel):
         lcp = self.createLCP()
         lcp.d[droppedlabel - 1] = 0  # subsidize this label
         # tabl.runlemke(verbose=True, lexstats=True, z0=gz0)
-        equilibrium = lemke.runlemke(lcp=lcp, silent=True)[1: lcp.n - 1]
+        equilibrium = lemke.runlemke(lcp=lcp)[1: lcp.n - 1]
         return tuple(equilibrium)
 
     def LH(self, LHstring):
@@ -279,7 +279,7 @@ def runtrace(self, xprior, yprior):
         Ay = self.A.negmatrix @ yprior
         xB = xprior @ self.B.negmatrix
         lcp.d = np.hstack((Ay, xB, [1, 1]))
-        equilibrium = lemke.runlemke(lcp=lcp, silent=True)[1: lcp.n - 1]
+        equilibrium = lemke.runlemke(lcp=lcp)[1: lcp.n - 1]
         return tuple(equilibrium)
 
     def tracing(self, trace):
diff --git a/src/lemke/lemke.py b/src/lemke/lemke.py
index cee3427..73ba08c 100644
--- a/src/lemke/lemke.py
+++ b/src/lemke/lemke.py
@@ -9,7 +9,6 @@
 # global defaults
 lcpfilename = "lcp"
 outfile = lcpfilename + ".out"
-filehandle = sys.stdout
 verbose = False
 silent = False
 z0 = False
@@ -45,10 +44,6 @@ def processArguments():
     return
 
 
-def printout(*s):
-    print(*s, file=filehandle)
-
-
 # LCP data M,q,d
 class lcp:
     # create LCP either with given n or from file
@@ -416,10 +411,6 @@ def __init__(self, enter, tableau):
         )
 
 
-def docupivot(leave, enter):  # leave, enter in VARS
-    return f"leaving: {leave.ljust(5)} entering: {enter}"
-
-
 def outsol(tableau):  # string giving solution, after createsol()
     # printout in columns to check complementarity
     n = tableau.n
@@ -468,8 +459,16 @@ def outstatistics(tableau):
     return stats
 
 
-def runlemke(*, lcp, verbose=False, lexstats=False, z0=False, silent=False):
-    global filehandle
+class LemkeCallback:
+    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
+
+
+class PrintingCallback(LemkeCallback):
     # z0: printout value of z0
     # flags.maxcount   = 0;
     # flags.bdocupivot = 1;
@@ -479,72 +478,117 @@ def runlemke(*, lcp, verbose=False, lexstats=False, z0=False, silent=False):
     # flags.binteract  = 0;
     # flags.blexstats  = 0;
 
+    def __init__(
+        self,
+        stream=sys.stdout,
+        verbose=False,
+        z0=False,
+        lexstats=False,
+    ):
+        self.stream = stream
+        self.verbose = verbose
+        self.z0 = z0
+        self.lexstats = lexstats
+
+    def printout(self, *args):
+        print(*args, file=self.stream)
+
+    def on_start(self, lcp, tableau):
+        self.printout(f"verbose={verbose} lcpfilename={lcpfilename} silent={silent} z0={z0}")
+        self.printout(lcp)
+        self.printout("==================================")
+
+        # if (flags.binitabl)
+        self.printout("After filltableau:")
+        self.printout(tableau)
+
+    def on_negcol(self, tableau):
+        # if (flags.binitabl)
+        if self.verbose:
+            self.printout("After negcol:")
+            self.printout(tableau)
+
+    def on_pivot_start(self, tableau, leave, enter):
+        if self.z0:  # printout progress of z0
+            z0_value = 0.0
+            if tableau.bascobas[0] < tableau.n:  # z0 is basic
+                z0_value = tableau.A[tableau.bascobas[0]][tableau.n + 1] / tableau.determinant
+            self.printout(f"pivot count = {tableau.pivotcount}, z0 = {z0_value}")
+
+        # if (flags.bdocupivot)
+        self.printout(f"leaving: {leave.ljust(5)} entering: {enter}")
+
+    def on_pivot_end(self, tableau):
+        if self.verbose:
+            self.printout(tableau)
+
+    def on_done(self, tableau):
+        if self.z0:
+            self.printout(f"pivot count = {tableau.pivotcount + 1}, z0 = 0.0")
+
+        # if (flags.binitabl)
+        self.printout("Final tableau:")
+        self.printout(tableau)
+
+        # if (flags.boutsol)
+        self.printout(outsol(tableau))
+
+        if self.lexstats:
+            # output statistics of minimum ratio test
+            self.printout(outstatistics(tableau))
+
+    def on_ray_termination(self, tableau, message):
+        self.printout(message)
+        self.printout(tableau)
+        self.printout("Current basis not an LCP solution:")
+        self.printout(outsol(tableau))
+
+
+def runlemke(*, lcp, callback=None):
+    callback = callback or LemkeCallback()
+
     try:
         tabl = tableau(lcp)
 
-        printout(f"verbose={verbose} lcpfilename={lcpfilename} silent={silent} z0={z0}")
-        printout(lcp)
-        printout("==================================")
-
-        if silent:
-            filehandle = open(outfile, "w")  # noqa: SIM115
         n = tabl.n
         tabl.pivotcount = 1
         # check if d is ok - TBC
         # if (flags.binitabl)
-        printout("After filltableau:")
-        printout(tabl)
+        callback.on_start(lcp=lcp, tableau=tabl)
 
         # z0 enters the basis to obtain lex-feasible solution
         enter = 0
         leave, z0leave = tabl.lexminvar(enter)
         # negate RHS
         tabl.negcol(n + 1)
-        # if (flags.binitabl)
-        if verbose:
-            printout("After negcol:")
-            printout(tabl)
+        callback.on_negcol(tableau=tabl)
 
         while True:  # main loop of complementary pivoting
             tabl.testtablvars()
-            if z0:  # printout progress of z0
-                if tabl.bascobas[0] < n:  # z0 is basic
-                    printout(
-                        "step,z0=",
-                        tabl.pivotcount, tabl.A[tabl.bascobas[0]][n + 1] / tabl.determinant
-                    )
-                else:
-                    printout("step,z0=", tabl.pivotcount, 0.0)
-            # if (flags.bdocupivot)
             tabl.assertbasic(leave, "docupivot")
             tabl.assertcobasic(enter, "docupivot")
-            printout(docupivot(tabl.vartoa(leave), tabl.vartoa(enter)))
+
+            callback.on_pivot_start(
+                tableau=tabl,
+                leave=tabl.vartoa(leave),
+                enter=tabl.vartoa(enter),
+            )
             tabl.pivot(leave, enter)
             if z0leave:
-                if z0:
-                    printout("step,z0=", tabl.pivotcount + 1, 0.0)
                 break
-            if verbose:
-                printout(tabl)
+
+            callback.on_pivot_end(tableau=tabl)
+
             enter = tabl.complement(leave)
             leave, z0leave = tabl.lexminvar(enter)
             tabl.pivotcount += 1
 
-        # if (flags.binitabl)
-        printout("Final tableau:")
-        printout(tabl)
-        # if (flags.boutsol)
         tabl.createsol()
-        printout(outsol(tabl))
-        if lexstats:
-            printout(outstatistics(tabl))
+        callback.on_done(tableau=tabl)
 
         return tabl.solution
     except RayTermination as e:
-        printout(str(e))
-        printout(e.tableau)
-        printout("Current basis not an LCP solution:")
-        printout(outsol(e.tableau))
+        callback.on_ray_termination(message=str(e), tableau=e.tableau)
         return None
 
 
@@ -553,7 +597,17 @@ def main():
 
     m = lcp(lcpfilename)
 
-    runlemke(lcp=m, verbose=verbose, z0=z0, silent=silent)
+    if silent:
+        with open(outfile, "w") as f:
+            runlemke(
+                lcp=m,
+                callback=PrintingCallback(stream=f, verbose=verbose, z0=z0),
+            )
+    else:
+        runlemke(
+            lcp=m,
+            callback=PrintingCallback(stream=sys.stdout, verbose=verbose, z0=z0),
+        )
 
 
 if __name__ == "__main__":

From a3c76db01ee994a7930c950fd49e080c5c1b8edc Mon Sep 17 00:00:00 2001
From: nataliemes <nataliemsbl@gmail.com>
Date: Fri, 24 Jul 2026 09:50:09 +0400
Subject: [PATCH 5/8] Update test_lemke.py

---
 tests/test_lemke.py | 11 ++++-------
 1 file changed, 4 insertions(+), 7 deletions(-)

diff --git a/tests/test_lemke.py b/tests/test_lemke.py
index f52b73d..1b66c0d 100644
--- a/tests/test_lemke.py
+++ b/tests/test_lemke.py
@@ -4,6 +4,7 @@
 from click.testing import CliRunner
 
 from lemke.lemke import (
+    RayTermination,
     lcp,
     main,
     tableau,
@@ -36,11 +37,9 @@ def test_lcp_invalid_file(tmp_path, content):
     file_path = tmp_path / "lcp"
     file_path.write_text(content)
 
-    with pytest.raises(SystemExit) as exc_info:
+    with pytest.raises(ValueError):
         lcp(str(file_path))
 
-    assert exc_info.value.code == 1
-
 
 # ---   TABLEAU INIT   ----------------------------------------------
 def test_tableau_dimensions():
@@ -90,7 +89,7 @@ def test_complement_pairs(i):
 
 def test_complement_z0_fails():
     t = tableau(lcp(2))
-    with pytest.raises(SystemExit):
+    with pytest.raises(RuntimeError):
         t.complement(0)
 
 
@@ -124,11 +123,9 @@ def test_lexminvar_without_positive_entry():
 
     enter = 0  # z0, cobasic in col 0
 
-    with pytest.raises(SystemExit) as exc_info:
+    with pytest.raises(RayTermination):
         t.lexminvar(enter)
 
-    assert exc_info.value.code == 1
-
 
 @pytest.mark.parametrize(
     "row0, row1",

From 4ab1cec2eac51a3d9449f61cba297107bec1cf24 Mon Sep 17 00:00:00 2001
From: nataliemes <nataliemsbl@gmail.com>
Date: Fri, 24 Jul 2026 16:07:48 +0400
Subject: [PATCH 6/8] CLI exit code 1 when no solution found

---
 src/lemke/lemke.py  |  5 ++++-
 tests/test_lemke.py | 12 +++++++++++-
 2 files changed, 15 insertions(+), 2 deletions(-)

diff --git a/src/lemke/lemke.py b/src/lemke/lemke.py
index db40b48..b72c08c 100644
--- a/src/lemke/lemke.py
+++ b/src/lemke/lemke.py
@@ -584,11 +584,14 @@ def main(verbose, z0, lcpfilename):
 
     m = lcp(lcpfilename)
 
-    runlemke(
+    result = runlemke(
         lcp=m,
         callback=PrintingCallback(stream=sys.stdout, verbose=verbose, z0=z0),
     )
 
+    if result is None:
+        sys.exit(1)
+
 
 if __name__ == "__main__":
     # m = lcp(3)
diff --git a/tests/test_lemke.py b/tests/test_lemke.py
index 1b66c0d..3961f89 100644
--- a/tests/test_lemke.py
+++ b/tests/test_lemke.py
@@ -151,7 +151,7 @@ def test_lexminvar_with_positive_entry(row0, row1):
     [],
     ["--verbose", "--z0"],
 ])
-def test_cli_runs_without_error(tmp_path, extra_args):
+def test_exit_code_0_on_success(tmp_path, extra_args):
     file_path = tmp_path / "test.lcp"
     file_path.write_text("n= 2\nM= 1 0 0 1\nq= 1 1\nd= 1 1\n")
 
@@ -161,6 +161,16 @@ def test_cli_runs_without_error(tmp_path, extra_args):
     assert result.exit_code == 0
 
 
+def test_exit_code_1_on_ray_termination(tmp_path):
+    file_path = tmp_path / "test.lcp"
+    file_path.write_text("n= 2\nM= -1 0 0 -1\nq= -1 -1\nd= 1 1\n")
+
+    runner = CliRunner()
+    result = runner.invoke(main, [str(file_path)])
+
+    assert result.exit_code == 1
+
+
 def test_cli_rejects_missing_file(tmp_path):
     missing_path = tmp_path / "missing"
 

From d7a2659ae7c353c509a99bb042a759afc67b546a Mon Sep 17 00:00:00 2001
From: nataliemes <nataliemsbl@gmail.com>
Date: Fri, 24 Jul 2026 16:10:12 +0400
Subject: [PATCH 7/8] Raise error in bimatrix if runlemke() finds no solution

---
 src/lemke/bimatrix.py | 14 ++++++++++++--
 1 file changed, 12 insertions(+), 2 deletions(-)

diff --git a/src/lemke/bimatrix.py b/src/lemke/bimatrix.py
index e7cdb20..cb518af 100644
--- a/src/lemke/bimatrix.py
+++ b/src/lemke/bimatrix.py
@@ -252,7 +252,12 @@ def runLH(self, droppedlabel):
         lcp = self.createLCP()
         lcp.d[droppedlabel - 1] = 0  # subsidize this label
         # tabl.runlemke(verbose=True, lexstats=True, z0=gz0)
-        equilibrium = lemke.runlemke(lcp=lcp)[1: lcp.n - 1]
+
+        result = lemke.runlemke(lcp=lcp)
+        if result is None:
+            raise RuntimeError("runlemke() failed to find a solution unexpectedly.")
+
+        equilibrium = result[1: lcp.n - 1]
         return tuple(equilibrium)
 
     def LH(self, LHstring):
@@ -279,7 +284,12 @@ def runtrace(self, xprior, yprior):
         Ay = self.A.negmatrix @ yprior
         xB = xprior @ self.B.negmatrix
         lcp.d = np.hstack((Ay, xB, [1, 1]))
-        equilibrium = lemke.runlemke(lcp=lcp)[1: lcp.n - 1]
+
+        result = lemke.runlemke(lcp=lcp)
+        if result is None:
+            raise RuntimeError("runlemke() failed to find a solution unexpectedly.")
+
+        equilibrium = result[1: lcp.n - 1]
         return tuple(equilibrium)
 
     def tracing(self, trace):

From 1533fc4f142802e4abfda1acd91a6e00ab8aef06 Mon Sep 17 00:00:00 2001
From: nataliemes <nataliemsbl@gmail.com>
Date: Fri, 31 Jul 2026 12:18:09 +0400
Subject: [PATCH 8/8] Update a test docstring

---
 tests/test_lcp.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_lcp.py b/tests/test_lcp.py
index a713d29..8b2d4e4 100644
--- a/tests/test_lcp.py
+++ b/tests/test_lcp.py
@@ -235,7 +235,7 @@ def test_with_lcp_conditions(test_case: LCPTestCase, subtests):
 def test_failure(test_case: LCPTestCase):
     """
     Test the Lemke solver on LCPs that terminate on a secondary ray
-    by verifying that it raises SystemExit with code 1.
+    by verifying that it returns None.
     """
     lcp_instance = test_case.factory()