Implementation — reading a minimal simplex

Two variables let you look at every vertex candidate

As Chapter 4 showed, the optimum normally appears at a vertex. In two dimensions, that means you can solve the LP directly with the following recipe.

1. List the boundary lines x = 0, y = 0, 2x + y = 16, x + 2y = 14.
2. Compute the intersection of every pair.
3. Keep only the intersections that satisfy every constraint.
4. Evaluate the objective at each surviving point and pick the maximum.

This method does not scale to general large-scale LPs, but it is perfect for an introduction.

Mapping the hand calculations to the code

The Python code you are about to read simply automates the work you already did by hand. Look at the correspondence first.

So the code is just Chapter 3 and Chapter 4 wrapped in a "pick two boundary lines" loop. The hard part is not the Python syntax but reading "how the manual procedure was laid out mechanically".

Why use a dict to represent each boundary line?

The code represents each boundary line as a dictionary: {"kind": "x0"}, {"kind": "y0"}, or {"kind": "line", "a": ..., "b": ..., "c": ...}. The reasons for using a dict instead of a class are:

• Only three kinds of boundary line exist: x = 0 (the y-axis), y = 0 (the x-axis), and a general line ax + by = c. A single string field kind is enough to tell them apart.
• Keep external dependencies at zero: by sticking to built-ins instead of dataclasses or custom classes, the code stays copy-and-runnable for the reader.
• Case-by-case dispatch is easy to read: "intersection with an axis" and "intersection of two general lines" use different formulas, so branching on kind ends up close to the equations themselves.

In a real codebase you would reach for @dataclass or typing.Literal; in a teaching example, making the case analysis visible is a feature, not a bug.

A minimal Python implementation

Without any external library, a minimal implementation that relies only on boundary intersections and a feasibility check looks like this.

EPS = 1e-9  # floating-point tolerance


def intersection(l1, l2):
    kind1 = l1["kind"]
    kind2 = l2["kind"]

    if kind1 == "x0" and kind2 == "y0":
        return (0.0, 0.0)
    if kind1 == "y0" and kind2 == "x0":
        return (0.0, 0.0)

    if kind1 == "x0" and kind2 == "line":
        return (0.0, l2["c"] / l2["b"])
    if kind1 == "line" and kind2 == "x0":
        return (0.0, l1["c"] / l1["b"])

    if kind1 == "y0" and kind2 == "line":
        return (l2["c"] / l2["a"], 0.0)
    if kind1 == "line" and kind2 == "y0":
        return (l1["c"] / l1["a"], 0.0)

    det = l1["a"] * l2["b"] - l2["a"] * l1["b"]
    if abs(det) < EPS:
        return None  # parallel lines: no unique intersection

    x = (l1["c"] * l2["b"] - l2["c"] * l1["b"]) / det
    y = (l1["a"] * l2["c"] - l2["a"] * l1["c"]) / det
    return (x, y)


def feasible(point, constraints):
    x, y = point
    if x < -EPS or y < -EPS:
        return False  # non-negativity x, y >= 0
    return all(a * x + b * y <= c + EPS for a, b, c in constraints)


def solve_lp_2d(cx, cy, constraints):
    if len(constraints) != 2:
        raise ValueError(
            "this 2D solver expects exactly 2 constraints; got "
            + str(len(constraints))
        )

    boundaries = [
        {"kind": "x0"},
        {"kind": "y0"},
        {"kind": "line", "a": constraints[0][0], "b": constraints[0][1], "c": constraints[0][2]},
        {"kind": "line", "a": constraints[1][0], "b": constraints[1][1], "c": constraints[1][2]},
    ]

    candidates = []
    for i in range(len(boundaries)):
        for j in range(i + 1, len(boundaries)):
            p = intersection(boundaries[i], boundaries[j])
            if p is None:
                continue
            if feasible(p, constraints):
                candidates.append((round(p[0], 9), round(p[1], 9)))

    vertices = sorted(set(candidates))
    if not vertices:
        raise ValueError("feasible region is empty")

    best = max(vertices, key=lambda p: cx * p[0] + cy * p[1])
    best_value = cx * best[0] + cy * best[1]
    return best, best_value, vertices

If the preconditions break and the feasible region is empty, the function raises ValueError directly. For a 2-variable teaching implementation, raising an explicit exception is easier to track during learning than silently handling the edge case.

Code reading notes

What does the threshold EPS = 1e-9 mean?

In floating-point arithmetic, a value that should be exactly 0 can come out around 1e-15. EPS = 1e-9 is the "differences smaller than this are considered numerical noise" tolerance.

• abs(det) < EPS: treat the two lines as parallel (no intersection) when the determinant is essentially 0; this also avoids a division-by-zero for 1 / 0.
• x < -EPS: allow a slightly negative numerical error against x ≥ 0 (otherwise a vertex that should sit exactly on the axis would be wrongly rejected).
• a * x + b * y <= c + EPS: the same logic applied to other constraints, so points right on the boundary are not thrown away by rounding.

The value 1e-9 is a common choice: large enough to absorb typical rounding error, small enough not to interfere with the actual problem scale. If your coefficients are on the order of 10⁶ or more, raising the tolerance to 1e-6 is a reasonable adjustment — match it to the scale of the numbers you handle.

Why skip the case det = 0 (parallel lines)?

When the determinant det = a₁b₂ − a₂b₁ is 0, the two boundary lines are parallel (or coincident). Parallel lines have no intersection (or infinitely many), so there is no useful vertex candidate to add. Returning None and reading if p is None: continue in the caller skips the case cleanly.

Why does intersection have so many branches?

There are only three kinds of boundary line (x0, y0, line), so the combinations are at most 3 × 3 = 9 in principle, or 6 with symmetry — and the code lists them explicitly. Because only four boundary lines exist in two variables, the brute force is called only 4C2 = 6 times in total. With n variables or many more constraints, this approach breaks down. Read it with the awareness that this is a simplification only allowed because we are in 2D.

Why deduplicate with sorted(set(candidates))?

The reason for rounding with round(p[0], 9) before passing through set is that the same vertex can show up several times with slightly different floating-point coordinates. For example, (0.0, 0.0) can be produced both by x0 ∩ y0 and by x0 ∩ line (when the y-intercept of the second constraint happens to be 0). Rounding to 9 decimal places collapses the practical differences, and set drops the duplicates.

Validating the number of constraints

Checking len(constraints) != 2 matters because this implementation assumes 4 boundary lines. If 3 or more constraints are passed, boundaries would silently ignore the extras and return a wrong answer; raising an exception at the entrance makes the mistake obvious.

Why is the non-negativity constraint handled separately?

Inside feasible, the non-negativity constraint x, y ≥ 0 is checked explicitly rather than included in the constraints array. That is not because non-negativity is theoretically different from any other constraint — it is the same kind of half-plane. The reason is that x = 0 and y = 0 are already in the boundary list, so the constraints array is reserved for the process constraints to avoid duplication. This is a readability decision, not a theoretical distinction.

Key takeaways from this chapter

• In two variables, enumerating every boundary intersection is still small enough to handle.
• After computing intersections, always use the feasibility check to drop points outside the region.
• In the end, you only need to compare the objective at the surviving vertices.