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"
VIOLETT = "#7c3aed"     # Niveaulinien und Zielfunktionsgradient

# 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_niveaulinien(achse, besucht: list[np.ndarray]) -> None:
    """Eine Niveaulinie der Zielfunktion durch JEDE besuchte Ecke.

    Die Niveaus werden aus dem Weg gerechnet, nicht gewaehlt: Jede Ecke des
    Simplexwegs liegt damit auf genau einer Linie, und der Abstand der Linien
    ist der Fortschritt eines Schritts. Damit zeigt das Bild nicht nur, WELCHEN
    Weg der Simplex nimmt, sondern warum - er schiebt die Linie so weit nach
    aussen, wie das Polyeder es zulaesst."""
    x = np.linspace(-0.9, 11.4, 100)
    for punkt in besucht:
        niveau = float(C @ punkt)
        achse.plot(x, (niveau - C[0] * x) / C[1], linestyle=":",
                   color=VIOLETT, linewidth=1.3, alpha=0.9, zorder=2)
        # Beschriftung dort, wo die Linie im Bild noch Platz hat. Die
        # Nulllinie laeuft nahe am unteren Rand - ihre Beschriftung gehoert
        # darueber, sonst wird sie abgeschnitten.
        unten = niveau < 500
        x_text = 1.7 if unten else 10.1
        y_linie = (niveau - C[0] * x_text) / C[1]
        achse.text(x_text, y_linie + (0.16 if unten else -0.28),
                   f"Z = {euro(niveau)}", fontsize=7.5, color=VIOLETT,
                   ha="center", va="bottom" if unten else "top", zorder=6)

    # Der Gradient: die Richtung, in die sich die Linien verschieben lassen.
    start = np.array([2.6, 2.2])
    ende = start + C / np.linalg.norm(C) * 2.6
    achse.annotate("", xy=tuple(ende), xytext=tuple(start),
                   arrowprops=dict(arrowstyle="-|>", color=VIOLETT,
                                   linewidth=2.2), zorder=6)
    achse.text(ende[0] + 0.15, ende[1], "$\\mathbf{c} = (150; 250)$",
               fontsize=8.5, color=VIOLETT, va="center", zorder=6)


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]
    zeichne_niveaulinien(achse, besucht)
    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)