operations_research/bilder_04/erzeuge_polyeder.py

#!/usr/bin/env python3

# erzeuge_polyeder.py
"""
Erzeugt die beiden Polyeder-Diagramme zum Fundament- und zum LP-Kapitel:

    bilder_04/kap_fundament_polyeder.svg   der zulaessige Bereich
    bilder_04/kap_lp_simplex_ecken.svg   der Eckenlauf des Simplex

Beide zeigen DASSELBE Modell - das Bot-Allokationsproblem, das sich durch die
ersten Kapitel zieht:

    max 150 x1 + 250 x2
    u.d.N.  2 x1 + 5 x2 <= 40   (vCPU)
            4 x1 + 6 x2 <= 60   (RAM)
                  x1    <=  8   (Marktlimit)
            x1, x2 >= 0

Warum das wichtig ist: Die Vorgaengerbilder zeigten ein voellig anderes LP
(2x1+3x2 <= 24, 3x1+2x2 <= 24, Optimum bei (4,8; 4,8)) - unmittelbar neben einer
Handrechnung, die mit dem Bot-Modell rechnet und (7,5; 5) herausbekommt. Ein
Leser musste annehmen, beides gehoere zusammen.

Die Ecken werden NICHT hingeschrieben, sondern nach demselben Verfahren
bestimmt wie in Visualisierung_Loesungsraum.py: Schnittpunkte je zweier
Begrenzungsgeraden, die zulaessig sind. Der Simplex-Weg im zweiten Bild ist der
Weg der Handrechnung im LP-Kapitel: (0,0) -> (0,8) -> (7,5; 5), also drei der
fuenf Ecken.

Aufruf (aus dem Repository-Wurzelverzeichnis):
    python3 bilder_04/erzeuge_polyeder.py

Benoetigt: numpy, matplotlib
"""

from __future__ import annotations

import itertools
import os

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np

# Reproduzierbare SVG-Ausgabe (siehe erzeuge_titelseite.py)
plt.rcParams["svg.hashsalt"] = "or-mit-python-v04"

BASIS = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
BILDER = os.path.join(BASIS, "bilder_04")

INDIGO, CYAN, AMBER, GRUEN = "#4338ca", "#0891b2", "#b45309", "#059669"

# Modell wortgleich zu Visualisierung_Loesungsraum.py
C = np.array([150.0, 250.0])
A = np.array([[2.0, 5.0], [4.0, 6.0], [1.0, 0.0]])
B = np.array([40.0, 60.0, 8.0])
NAMEN = ["vCPU: $2x_1+5x_2 \\leq 40$", "RAM: $4x_1+6x_2 \\leq 60$",
         "Marktlimit: $x_1 \\leq 8$"]
FARBEN = [INDIGO, CYAN, AMBER]

# Der Weg der Handrechnung im LP-Kapitel: Start im Ursprung, dann x2 aufnehmen
# (Pivotspalte -250), dann x1 (Pivotspalte -50).
SIMPLEXWEG = [(0.0, 0.0), (0.0, 8.0), (7.5, 5.0)]


def ist_zulaessig(punkt: np.ndarray, tol: float = 1e-7) -> bool:
    return bool(np.all(A @ punkt <= B + tol) and np.all(punkt >= -tol))


def berechne_ecken() -> np.ndarray:
    """Ecken = zulaessige Schnittpunkte je zweier Begrenzungsgeraden.

    Wortgleich zu Visualisierung_Loesungsraum.py - die Koordinaten im Bild
    stammen also aus derselben Rechnung wie die im Buch abgedruckte Tabelle,
    nicht aus einer zweiten Quelle.
    """
    geraden = [(A[i], B[i]) for i in range(len(B))]
    geraden.append((np.array([1.0, 0.0]), 0.0))     # x1 = 0
    geraden.append((np.array([0.0, 1.0]), 0.0))     # x2 = 0

    ecken: list[np.ndarray] = []
    for (n1, d1), (n2, d2) in itertools.combinations(geraden, 2):
        matrix = np.array([n1, n2])
        if abs(np.linalg.det(matrix)) < 1e-9:
            continue
        punkt = np.linalg.solve(matrix, np.array([d1, d2]))
        if ist_zulaessig(punkt) and not any(np.allclose(punkt, e) for e in ecken):
            ecken.append(punkt)
    return np.array(ecken)


# Beschriftungsposition je Ecke, in Punkten relativ zum Eckpunkt. Von Hand
# gesetzt und nicht gerechnet: (7,5; 5) und (8; 4,67) liegen eine halbe Einheit
# auseinander, und jede automatische Regel legt ihre Beschriftungen uebereinander.
# Das ist Layout, keine Aussage - die Koordinaten selbst bleiben gerechnet.
VERSATZ = {
    (0.0, 0.0):  (14, 10),
    (0.0, 8.0):  (14, -2),
    (8.0, 0.0):  (-10, 14),
    (7.5, 5.0):  (-12, 16),
    (8.0, 4.67): (16, -4),
}


# Im Simplex-Bild laeuft der Pfad durch zwei dieser Positionen - dort andere.
VERSATZ_SIMPLEX = {
    (0.0, 8.0): (16, 12),
    (7.5, 5.0): (14, 16),
}


def versatz(punkt: np.ndarray, mit_pfad: bool = False) -> tuple:
    schluessel = (round(float(punkt[0]), 2), round(float(punkt[1]), 2))
    if mit_pfad and schluessel in VERSATZ_SIMPLEX:
        return VERSATZ_SIMPLEX[schluessel]
    return VERSATZ.get(schluessel, (12, 8))


def punkt_text(punkt: np.ndarray) -> str:
    """Koordinaten deutsch und auf zwei Stellen - '(8; 4,67)' statt '(8; 4.66667)'."""
    teile = []
    for wert in punkt:
        text = f"{wert:.2f}".rstrip("0").rstrip(".")
        teile.append(text.replace(".", ","))
    return f"({teile[0]}; {teile[1]})"


def euro(wert: float) -> str:
    return f"{wert:,.0f} €".replace(",", ".")


def sortiere_umlaufend(ecken: np.ndarray) -> np.ndarray:
    """Ecken gegen den Uhrzeigersinn ordnen, damit das Polygon nicht verknotet."""
    mitte = ecken.mean(axis=0)
    winkel = np.arctan2(ecken[:, 1] - mitte[1], ecken[:, 0] - mitte[0])
    return ecken[np.argsort(winkel)]


def zeichne_rahmen(achse, ecken: np.ndarray) -> None:
    """Restriktionsgeraden, zulaessiger Bereich und Achsenbeschriftung."""
    rand = sortiere_umlaufend(ecken)
    achse.fill(rand[:, 0], rand[:, 1], color=INDIGO, alpha=0.10, zorder=1)
    achse.plot(np.append(rand[:, 0], rand[0, 0]),
               np.append(rand[:, 1], rand[0, 1]),
               color=INDIGO, linewidth=1.6, zorder=2)

    x = np.linspace(0, 11, 200)
    for zeile, grenze, name, farbe in zip(A, B, NAMEN, FARBEN):
        if abs(zeile[1]) < 1e-9:                     # senkrechte Gerade x1 = c
            achse.axvline(grenze / zeile[0], color=farbe, linewidth=1.3,
                          linestyle="--", alpha=0.8, label=name)
        else:
            achse.plot(x, (grenze - zeile[0] * x) / zeile[1], color=farbe,
                       linewidth=1.3, linestyle="--", alpha=0.8, label=name)

    achse.set_xlim(-0.9, 11.4)
    achse.set_ylim(-0.9, 11)
    achse.set_xlabel("$x_1$ — Bots vom Typ A")
    achse.set_ylabel("$x_2$ — Bots vom Typ B")
    achse.grid(linestyle=":", alpha=0.45)
    for seite in ("top", "right"):
        achse.spines[seite].set_visible(False)


def schreibe(figur, name: str) -> None:
    os.makedirs(BILDER, exist_ok=True)
    for endung in ("svg",):
        pfad = os.path.join(BILDER, f"{name}.{endung}")
        figur.savefig(pfad, format=endung, dpi=160,
                      metadata={"Date": None} if endung == "svg" else None)
        print(f"geschrieben: {pfad}")
    plt.close(figur)


def zeichne_polyeder(ecken: np.ndarray, bestes: np.ndarray) -> None:
    """Fundament: alle Ecken mit ihrem Zielwert - der Fundamentalsatz zum Anfassen."""
    figur, achse = plt.subplots(figsize=(7.2, 5.2))
    zeichne_rahmen(achse, ecken)

    for ecke in ecken:
        wert = float(C @ ecke)
        optimal = np.allclose(ecke, bestes)
        achse.scatter(*ecke, s=120 if optimal else 55, zorder=4,
                      color=GRUEN if optimal else INDIGO)
        beschriftung = f"{punkt_text(ecke)}\n{euro(wert)}"
        if optimal:
            beschriftung += "\n← Optimum"
        weg = versatz(ecke)
        achse.annotate(beschriftung, ecke, textcoords="offset points",
                       xytext=weg, fontsize=8.5,
                       ha="left" if weg[0] >= 0 else "right",
                       va="bottom" if weg[1] >= 0 else "top",
                       color=GRUEN if optimal else "#334155",
                       fontweight="bold" if optimal else "normal")
    achse.set_title("Der zulässige Bereich und seine Ecken\n"
                    "$\\max\\ 150x_1 + 250x_2$", fontsize=11)
    achse.legend(fontsize=8, loc="upper right", frameon=False)
    figur.tight_layout()
    schreibe(figur, "kap_fundament_polyeder")


def zeichne_simplexweg(ecken: np.ndarray, bestes: np.ndarray) -> None:
    """LP-Kapitel: derselbe Bereich, aber der Weg der Handrechnung."""
    figur, achse = plt.subplots(figsize=(7.2, 5.2))
    zeichne_rahmen(achse, ecken)

    besucht = [np.array(p) for p in SIMPLEXWEG]
    unbesucht = [e for e in ecken
                 if not any(np.allclose(e, p) for p in besucht)]
    for ecke in unbesucht:
        achse.scatter(*ecke, s=55, color="#94a3b8", zorder=3)
        weg = versatz(ecke, mit_pfad=True)
        achse.annotate(punkt_text(ecke), ecke, textcoords="offset points",
                       xytext=weg, fontsize=8, color="#64748b",
                       ha="left" if weg[0] >= 0 else "right",
                       va="bottom" if weg[1] >= 0 else "top")

    weg = np.array(besucht)
    achse.plot(weg[:, 0], weg[:, 1], color=AMBER, linewidth=2.4, zorder=4)
    for schritt, punkt in enumerate(besucht):
        wert = float(C @ punkt)
        letzter = schritt == len(besucht) - 1
        achse.scatter(*punkt, s=140 if letzter else 90, zorder=5,
                      color=GRUEN if letzter else AMBER)
        weg = versatz(punkt, mit_pfad=True)
        achse.annotate(f"{schritt}. {punkt_text(punkt)}\nZ = {euro(wert)}",
                       punkt, textcoords="offset points", xytext=weg,
                       fontsize=8.5, fontweight="bold",
                       ha="left" if weg[0] >= 0 else "right",
                       va="bottom" if weg[1] >= 0 else "top",
                       color=GRUEN if letzter else "#334155")

    achse.set_title(f"Der Simplex besucht {len(besucht)} von {len(ecken)} Ecken\n"
                    "und keine davon zweimal", fontsize=11)
    achse.legend(fontsize=8, loc="upper right", frameon=False)
    figur.tight_layout()
    schreibe(figur, "kap_lp_simplex_ecken")


if __name__ == "__main__":
    ecken = berechne_ecken()
    werte = ecken @ C
    bestes = ecken[int(np.argmax(werte))]
    print(f"{len(ecken)} Ecken gefunden:")
    for ecke, wert in sorted(zip(ecken.tolist(), werte.tolist()),
                             key=lambda p: -p[1]):
        print(f"  ({ecke[0]:6.2f}; {ecke[1]:5.2f})  Z = {wert:9,.2f} EUR")
    print(f"Optimum: ({bestes[0]:g}; {bestes[1]:g})")
    for punkt in SIMPLEXWEG:
        if not any(np.allclose(np.array(punkt), e) for e in ecken):
            raise SystemExit(f"FEHLER: {punkt} ist keine Ecke des Polyeders.")
    zeichne_polyeder(ecken, bestes)
    zeichne_simplexweg(ecken, bestes)