operations_research/bilder_04/erzeuge_kondition.py

#!/usr/bin/env python3

# erzeuge_kondition.py
"""
Erzeugt das Konditionszahl-Diagramm zum Kapitel Fundament:

    bilder_04/kap_fundament_kondition.svg   statisch, fuer das PDF
    bilder_04/kap_fundament_kondition.png   Rasterfassung

Zwei Aussagen, beide gerechnet und nicht behauptet:

  Links  - kappa(A) ist keine grobe Faustregel, sondern eine SCHARFE obere
           Schranke. Die gemessene groesste Fehlerverstaerkung liegt ueber
           zwoelf Groessenordnungen hinweg praktisch auf der Geraden y = kappa;
           der Median darunter. Wer also mit kappa = 1e8 rechnet, verliert im
           schlimmsten Fall genau acht Stellen - nicht ungefaehr, sondern
           nachweisbar.
  Rechts - Was die Ruiz-Equilibrierung an einem Modell ausrichtet, in dem
           Euro-Betraege und Tonnen-Angaben in derselben Matrix stehen.

Verfahren und Instanz sind wortgleich zu Skalierung_Kondition.py, und
gerechnet wird hier erneut - so koennen Diagramm und Buchtext nicht
auseinanderlaufen.

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

Benoetigt: numpy, matplotlib
"""

from __future__ import annotations

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"

STOERUNG = 1e-10          # relative Stoerung der rechten Seite

GUT = np.array([[2.0, 1.0], [1.0, 3.0]])
SCHLECHT = np.array([[1.0, 1.0], [1.0, 1.0 + 1e-8]])


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

def fehlerverstaerkung(A: np.ndarray, rng: np.random.Generator,
                       versuche: int = 200) -> tuple[float, float]:
    """Stoert b relativ um STOERUNG und misst den Faktor, um den der
    Loesungsfehler groesser ausfaellt. Liefert (Median, Maximum)."""
    x_wahr = np.ones(A.shape[1])
    b = A @ x_wahr
    faktoren = []
    for _ in range(versuche):
        richtung = rng.normal(size=b.size)
        richtung /= np.linalg.norm(richtung)
        b_gestoert = b + STOERUNG * np.linalg.norm(b) * richtung
        x_gestoert = np.linalg.solve(A, b_gestoert)
        rel_x = np.linalg.norm(x_gestoert - x_wahr) / np.linalg.norm(x_wahr)
        faktoren.append(rel_x / STOERUNG)
    return float(np.median(faktoren)), float(np.max(faktoren))


def ruiz_equilibrierung(A: np.ndarray, durchlaeufe: int = 20) -> np.ndarray:
    """Zeilen- und Spaltenskalierung, bis alle Maxima nahe 1 liegen."""
    A_s = A.astype(float).copy()
    for _ in range(durchlaeufe):
        zeilen_max = np.abs(A_s).max(axis=1)
        zeilen_max[zeilen_max == 0] = 1.0
        A_s = (1.0 / np.sqrt(zeilen_max))[:, None] * A_s
        spalten_max = np.abs(A_s).max(axis=0)
        spalten_max[spalten_max == 0] = 1.0
        A_s = A_s * (1.0 / np.sqrt(spalten_max))[None, :]
    return A_s


def mischmodell() -> np.ndarray:
    """Das Modell mit den vier unvertraeglichen Einheiten.

    Skalierung_Kondition.py zieht diese Matrix aus DEMSELBEN Generator, aus dem
    zuvor schon die Stoerungen von Experiment 1 kamen. Damit hier bitgleich
    dieselbe Matrix entsteht - und das Diagramm dieselbe Konditionszahl zeigt
    wie der abgedruckte Programmlauf -, wird die Ziehungsreihenfolge exakt
    nachgespielt: erst 2 x 200 Stoerungen, dann die Grundmatrix.
    """
    rng = np.random.default_rng(42)
    for A in (GUT, SCHLECHT):
        fehlerverstaerkung(A, rng)
    grundmatrix = rng.uniform(0.5, 2.0, size=(4, 4))
    einheiten = np.array([1e7, 1e-3, 1e5, 1e-2])
    return grundmatrix * einheiten[:, None]


def messreihe() -> tuple[list[float], list[float], list[float]]:
    """Eine Familie fast paralleler Zeilen: A = [[1, 1], [1, 1 + eps]].

    Mit fallendem eps waechst kappa(A) um zwoelf Groessenordnungen, ohne dass
    sich sonst irgendetwas am Modell aendert. Genau das macht die Kurve
    aussagekraeftig: Der einzige Unterschied ist die Kondition.
    """
    # Eigener Generator, damit die Messreihe den nachgespielten Strom in
    # mischmodell() nicht stoert.
    rng = np.random.default_rng(7)
    kappas, mediane, maxima = [], [], []
    for exponent in np.arange(-1.0, -12.5, -0.5):
        eps = 10.0 ** exponent
        A = np.array([[1.0, 1.0], [1.0, 1.0 + eps]])
        median, maximum = fehlerverstaerkung(A, rng)
        kappas.append(float(np.linalg.cond(A)))
        mediane.append(median)
        maxima.append(maximum)
    return kappas, mediane, maxima


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

def zeichne(kappas, mediane, maxima, kappa_vorher, kappa_nachher) -> None:
    figur, (links, rechts) = plt.subplots(
        1, 2, figsize=(10.4, 4.3), gridspec_kw={"width_ratios": [1.75, 1]})

    # --- links: Verstaerkung gegen kappa ---
    # Die Gerade reicht bewusst ueber die Messpunkte hinaus: Innerhalb des
    # Messbereichs liegt sie exakt unter der Messkurve - und genau das ist die
    # Aussage. Sichtbar bleibt sie dadurch an den Enden.
    grenze = np.array([min(kappas) / 30, max(kappas) * 30])
    links.plot(grenze, grenze, color="#94a3b8", linewidth=3.0, linestyle="-",
               alpha=0.55, zorder=1,
               label="theoretische Schranke $y=\\kappa(A)$")
    links.plot(kappas, maxima, "o-", color=INDIGO, markersize=4.5, linewidth=1.4,
               zorder=3, label="gemessen: größte Verstärkung")
    links.plot(kappas, mediane, "s-", color=CYAN, markersize=4.0, linewidth=1.4,
               zorder=2, label="gemessen: Median")

    marken_rng = np.random.default_rng(7)
    for A, name, versatz in ((GUT, "gut konditioniert", (14, -4)),
                             (SCHLECHT, "fast parallele Zeilen", (-12, 14))):
        _, maximum = fehlerverstaerkung(A, marken_rng)
        k = float(np.linalg.cond(A))
        links.scatter([k], [maximum], s=95, facecolors="none", edgecolors=AMBER,
                      linewidths=2.0, zorder=4)
        links.annotate(name, (k, maximum), textcoords="offset points",
                       xytext=versatz, fontsize=8, color=AMBER,
                       ha="right" if versatz[0] < 0 else "left")

    links.set_xscale("log")
    links.set_yscale("log")
    links.set_xlabel("Konditionszahl $\\kappa(A)$")
    links.set_ylabel("Fehlerverstärkung (Faktor)")
    links.set_title("$\\kappa(A)$ ist eine scharfe Schranke, keine Faustregel",
                    fontsize=10.5)
    links.grid(which="both", linestyle=":", alpha=0.45)
    links.legend(fontsize=8, loc="upper left", frameon=False)
    links.annotate("Die Messung liegt auf der Schranke:\n"
                   "$\\kappa$ wird im schlimmsten Fall\nvollständig ausgeschöpft.",
                   xy=(0.97, 0.06), xycoords="axes fraction", fontsize=8,
                   ha="right", va="bottom", color="#475569")
    for rand in ("top", "right"):
        links.spines[rand].set_visible(False)

    # --- rechts: was Ruiz ausrichtet ---
    balken = rechts.bar(["vorher", "nach Ruiz"], [kappa_vorher, kappa_nachher],
                        color=[AMBER, GRUEN], width=0.55)
    rechts.set_yscale("log")
    rechts.set_ylim(1, kappa_vorher * 40)
    for stab, wert in zip(balken, (kappa_vorher, kappa_nachher)):
        rechts.annotate(f"{wert:.2e}", (stab.get_x() + stab.get_width() / 2, wert),
                        textcoords="offset points", xytext=(0, 6),
                        ha="center", fontsize=9, fontweight="bold")
    stellen = int(round(np.log10(kappa_vorher / kappa_nachher)))
    rechts.set_title(f"Ruiz-Equilibrierung:\n{stellen} Größenordnungen besser",
                     fontsize=10.5)
    rechts.set_ylabel("$\\kappa(A)$")
    rechts.grid(axis="y", which="both", linestyle=":", alpha=0.45)
    for rand in ("top", "right"):
        rechts.spines[rand].set_visible(False)

    figur.tight_layout()
    os.makedirs(BILDER, exist_ok=True)
    for endung in ("svg", "png"):
        pfad = os.path.join(BILDER, f"kap_fundament_kondition.{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__":
    kappas, mediane, maxima = messreihe()
    modell = mischmodell()
    kappa_vorher = float(np.linalg.cond(modell))
    kappa_nachher = float(np.linalg.cond(ruiz_equilibrierung(modell)))
    print(f"Messreihe: {len(kappas)} Matrizen, kappa von {min(kappas):.2e} "
          f"bis {max(kappas):.2e}")
    print(f"Ruiz: kappa {kappa_vorher:.2e} -> {kappa_nachher:.2e}")
    zeichne(kappas, mediane, maxima, kappa_vorher, kappa_nachher)