operations_research/bilder_04/erzeuge_runden.py

#!/usr/bin/env python3

# erzeuge_runden.py
"""
Erzeugt das Diagramm "Warum Runden scheitert" zum MILP-Kapitel:

    bilder_04/kap_milp_runden.svg   statisch, fuer das PDF

Zwei Bilder nebeneinander, beide aus dem Kapitel selbst:

  Links  - die Handrechnung des Kapitels, gezeichnet:
               max x1 + x2   u.d.N.  2 x1 + 2 x2 <= 3,  x ganzzahlig
           Der Witz ist die Entartung: Die Zielfunktion laeuft PARALLEL zur
           Restriktion, also ist die ganze Kante optimal, nicht eine Ecke.
           Aufrunden auf (1,1) fuehrt aus dem zulaessigen Bereich heraus,
           Abrunden auf (0,0) auf den Wert null. Beides sieht man sofort.
  Rechts - die Messung aus Runden_Gegenbeispiel.py, mit derselben Saat
           nachgerechnet: Wie oft ist die aufgerundete Loesung unzulaessig,
           und was kostet Abrunden - mit und ohne gieriges Auffuellen?

Warum ueberhaupt neu: Das Vorgaengerbild zeigte ein LP, das im Buch nicht
vorkommt (LP-Optimum bei (1,8; 3,4), Z = 42,6), ohne die Zielfunktion zu nennen
- unmittelbar ueber der Handrechnung, die etwas voellig anderes rechnet.

Aufruf (aus dem Repository-Wurzelverzeichnis, dauert rund eine Minute):
    python3 bilder_04/erzeuge_runden.py

Benoetigt: numpy, scipy, matplotlib
"""

from __future__ import annotations

import os

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import linprog

# 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, ROT = (
    "#4338ca", "#0891b2", "#b45309", "#059669", "#be123c")

GROESSEN = [(5, 2), (10, 3), (20, 5), (40, 8)]
INSTANZEN = 200
SAAT = 2026


# --- Messung: wortgleich zu Runden_Gegenbeispiel.py -------------------------

def erzeuge_instanz(n, m, rng):
    c = rng.integers(3, 20, size=n).astype(float)
    A = rng.integers(1, 9, size=(m, n)).astype(float)
    b = (A.sum(axis=1) * rng.uniform(0.25, 0.45)).round()
    return c, A, b


def loese_lp(c, A, b, ganzzahlig=False):
    n = len(c)
    ergebnis = linprog(
        c=-c, A_ub=A, b_ub=b, bounds=[(0, None)] * n,
        integrality=np.ones(n) if ganzzahlig else None,
        method="highs")
    return (-ergebnis.fun, ergebnis.x) if ergebnis.success else (None, None)


def abrunden_und_reparieren(x_lp, c, A, b):
    x = np.floor(x_lp + 1e-9)
    verbessert = True
    while verbessert:
        verbessert = False
        for j in np.argsort(-c):
            kandidat = x.copy()
            kandidat[j] += 1
            if np.all(A @ kandidat <= b + 1e-9):
                x = kandidat
                verbessert = True
                break
    return c @ x, x


def messreihe():
    """Dieselbe Schleife wie im Buchprogramm, mit derselben Saat.

    Die Saat sitzt VOR der Schleife ueber die Groessen, nicht darin - so wie
    dort. Nur deshalb kommen hier dieselben Zahlen heraus wie in der
    abgedruckten Tabelle.
    """
    rng = np.random.default_rng(SAAT)
    zeilen = []
    for n, m in GROESSEN:
        unzulaessig = 0
        verluste_ab, verluste_gierig = [], []
        gezaehlt = 0
        for _ in range(INSTANZEN):
            c, A, b = erzeuge_instanz(n, m, rng)
            z_lp, x_lp = loese_lp(c, A, b, ganzzahlig=False)
            z_ip, _ = loese_lp(c, A, b, ganzzahlig=True)
            if z_lp is None or z_ip is None or z_ip <= 0:
                continue
            gezaehlt += 1
            x_auf = np.ceil(x_lp - 1e-9)
            if np.any(A @ x_auf > b + 1e-9):
                unzulaessig += 1
            x_ab = np.floor(x_lp + 1e-9)
            verluste_ab.append(1.0 - (c @ x_ab) / z_ip)
            z_gierig, _ = abrunden_und_reparieren(x_lp, c, A, b)
            verluste_gierig.append(1.0 - z_gierig / z_ip)
        zeilen.append({
            "name": f"{n}×{m}",
            "unzulaessig": unzulaessig / gezaehlt,
            "ab": float(np.mean(verluste_ab)),
            "gierig": float(np.mean(verluste_gierig)),
        })
    return zeilen


# --- Zeichnen ---------------------------------------------------------------

def zeichne_handrechnung(achse) -> None:
    """Die Handrechnung des Kapitels: max x1+x2 u.d.N. 2x1+2x2 <= 3."""
    ecken = np.array([[0.0, 0.0], [1.5, 0.0], [0.0, 1.5]])
    achse.fill(ecken[:, 0], ecken[:, 1], color=CYAN, alpha=0.13, zorder=1)
    achse.plot([1.5, 0.0], [0.0, 1.5], color=CYAN, linewidth=2.0, zorder=3,
               label="$2x_1+2x_2 = 3$ — zugleich das LP-Optimum")

    # Ganzzahlige Punkte im Fenster, zulaessig und unzulaessig getrennt
    for x1 in range(0, 3):
        for x2 in range(0, 3):
            zulaessig = 2 * x1 + 2 * x2 <= 3
            achse.scatter(x1, x2, s=48, zorder=4,
                          color=INDIGO if zulaessig else "#cbd5e1")

    achse.scatter([0.75], [0.75], marker="*", s=260, color=ROT, zorder=6)
    achse.annotate("LP-Lösung $(0{,}75;\\ 0{,}75)$\n$Z_{LP} = 1{,}5$",
                   (0.75, 0.75), textcoords="offset points", xytext=(-12, -6),
                   ha="right", va="top", fontsize=9, color=ROT,
                   fontweight="bold")

    achse.scatter([1], [1], s=150, facecolors="none", edgecolors=AMBER,
                  linewidths=2.4, zorder=6)
    achse.annotate("aufgerundet $(1;1)$\n$2+2=4>3$ — unzulässig",
                   (1, 1), textcoords="offset points", xytext=(12, -6),
                   fontsize=8.5, color=AMBER, va="top")

    achse.scatter([0], [0], s=150, facecolors="none", edgecolors="#64748b",
                  linewidths=2.2, zorder=6)
    achse.annotate("abgerundet $(0;0)$\n$Z = 0$ — alles verschenkt",
                   (0, 0), textcoords="offset points", xytext=(14, 8),
                   fontsize=8.5, color="#475569")

    for punkt in ((1, 0), (0, 1)):
        achse.scatter(*punkt, s=140, color=GRUEN, zorder=7)
    achse.annotate("wahres Optimum $Z_{IP}=1$", (1, 0),
                   textcoords="offset points", xytext=(12, -4),
                   va="top", fontsize=8.5, color=GRUEN, fontweight="bold")

    achse.set_xlim(-0.35, 2.35)
    achse.set_ylim(-0.35, 2.35)
    achse.set_xlabel("$x_1$")
    achse.set_ylabel("$x_2$")
    achse.set_title("Die Handrechnung als Bild:\nrunden geht in beide Richtungen schief",
                    fontsize=10.5)
    achse.grid(linestyle=":", alpha=0.45)
    achse.legend(fontsize=8, loc="upper right", frameon=False)
    for seite in ("top", "right"):
        achse.spines[seite].set_visible(False)


def zeichne_messung(achse, zeilen) -> None:
    namen = [z["name"] for z in zeilen]
    stellen = np.arange(len(zeilen))
    breite = 0.26

    achse.bar(stellen - breite, [z["unzulaessig"] * 100 for z in zeilen],
              breite, color=AMBER, label="Aufrunden unzulässig")
    achse.bar(stellen, [z["ab"] * 100 for z in zeilen],
              breite, color=ROT, label="Abrunden: mittlerer Verlust")
    achse.bar(stellen + breite, [z["gierig"] * 100 for z in zeilen],
              breite, color=GRUEN, label="abrunden + gierig auffüllen")

    for stelle, zeile in zip(stellen, zeilen):
        for versatz, wert in ((-breite, zeile["unzulaessig"]),
                              (0.0, zeile["ab"]), (breite, zeile["gierig"])):
            achse.annotate(f"{wert * 100:.0f}", (stelle + versatz, wert * 100),
                           textcoords="offset points", xytext=(0, 3),
                           ha="center", fontsize=7.5, color="#334155")

    achse.set_xticks(stellen)
    achse.set_xticklabels(namen)
    achse.set_xlabel("Instanzgröße (Variablen × Nebenbedingungen)")
    achse.set_ylabel("Anteil bzw. Verlust (%)")
    achse.set_ylim(0, 118)
    achse.set_title(f"Gemessen an je {INSTANZEN} Zufallsinstanzen:\n"
                    "Aufrunden wird mit der Größe hoffnungslos", fontsize=10.5)
    achse.grid(axis="y", linestyle=":", alpha=0.45)
    achse.legend(fontsize=8, loc="upper center", frameon=False, ncol=1)
    for seite in ("top", "right"):
        achse.spines[seite].set_visible(False)


def zeichne(zeilen) -> None:
    figur, (links, rechts) = plt.subplots(1, 2, figsize=(11.0, 4.8))
    zeichne_handrechnung(links)
    zeichne_messung(rechts, zeilen)
    figur.tight_layout()

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


if __name__ == "__main__":
    zeilen = messreihe()
    print(f"{'n x m':>8} | {'Aufrunden unzul.':>17} | {'Abrunden':>10} | {'gierig':>8}")
    print("-" * 54)
    for zeile in zeilen:
        print(f"{zeile['name']:>8} | {zeile['unzulaessig'] * 100:>16.1f} % | "
              f"{zeile['ab'] * 100:>9.1f} % | {zeile['gierig'] * 100:>7.1f} %")
    zeichne(zeilen)