operations_research/bilder_04/erzeuge_predict_vs_optimize.py

#!/usr/bin/env python3

# erzeuge_predict_vs_optimize.py
"""
Erzeugt das Diagramm zum Predict-then-Optimize-Kapitel:

    bilder_04/kap_prognose_entscheidungsguete.svg

Vier Verfahren, zwei Massstaebe: links der Prognosefehler (MSE), rechts die
Kosten je Tag. Die Balken sind absichtlich nebeneinander gestellt, denn der
Befund des Kapitels ist genau ihr Auseinanderlaufen: Das Verfahren mit dem
BESTEN MSE hat die HOECHSTEN Kosten.

Ohne Formel gesagt: Der Prognostiker optimiert seinen Fehler, der Planer traegt
die Kosten - und beide messen etwas anderes. Wer die Prognose verbessert, ohne
zu fragen, wofuer sie gebraucht wird, verbessert die falsche Zahl.

Alle Werte kommen aus Predict_then_Optimize.py: dieselben Daten, dieselbe Saat,
dieselben vier Verfahren. pruefe_befund() stellt sicher, dass die Aussage des
Kapitels im Bild auch wirklich zu sehen ist - sonst waere das Bild ein
Widerspruch zum Text daneben.

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

Benoetigt: numpy, scipy, scikit-learn, matplotlib
"""

from __future__ import annotations

import os
import sys

import numpy as np

HIER = os.path.dirname(os.path.abspath(__file__))
BASIS = os.path.dirname(HIER)
sys.path.insert(0, HIER)
sys.path.insert(0, os.path.join(
    BASIS, "Operations_Research_mit_Python_Version_04_Programme"))

from stil_04 import FARBEN, speichere                          # noqa: E402
import matplotlib.pyplot as plt                                # noqa: E402
from scipy.stats import norm                                   # noqa: E402
from sklearn.linear_model import LinearRegression, QuantileRegressor  # noqa: E402
from Predict_then_Optimize import (KRITISCHES_VERHAELTNIS, TAGE,  # noqa: E402
                                   TRAINING, erzeuge_daten, tageskosten)

KURZ = ["Punkt-\nprognose", "+ Zuschlag aus der\nNormalverteilung",
        "+ Zuschlag auf\nKosten trainiert", "Quantil-\nregression"]


def rechne() -> list[tuple[str, float, float]]:
    """Dieselben vier Verfahren wie im Buchprogramm."""
    merkmale, nachfrage = erzeuge_daten()[:2]
    # Dieselbe Aufteilung wie im Buchprogramm - eine selbst gewaehlte Grenze
    # ergaebe andere Zahlen als die abgedruckte Tabelle.
    lernen = slice(0, TRAINING)
    pruefen = slice(TRAINING, TAGE)

    kq = LinearRegression().fit(merkmale[lernen], nachfrage[lernen])
    punkt = kq.predict(merkmale[pruefen])
    rest = nachfrage[lernen] - kq.predict(merkmale[lernen])
    pauschal = norm.ppf(KRITISCHES_VERHAELTNIS) * rest.std()

    kandidaten = np.linspace(-10.0, 30.0, 401)
    trainingsprognose = kq.predict(merkmale[lernen])
    kostenzuschlag = float(kandidaten[np.argmin(
        [tageskosten(trainingsprognose + z, nachfrage[lernen])
         for z in kandidaten])])

    quantil = QuantileRegressor(quantile=KRITISCHES_VERHAELTNIS, alpha=0.0,
                                solver="highs").fit(merkmale[lernen],
                                                    nachfrage[lernen])
    quantilprognose = quantil.predict(merkmale[pruefen])

    verfahren = [
        (KURZ[0], punkt, punkt),
        (KURZ[1], punkt, punkt + pauschal),
        (KURZ[2], punkt, punkt + kostenzuschlag),
        (KURZ[3], quantilprognose, quantilprognose),
    ]
    aus = []
    for name, prognose, bestellung in verfahren:
        mse = float(((prognose - nachfrage[pruefen]) ** 2).mean())
        aus.append((name, mse, tageskosten(bestellung, nachfrage[pruefen])))
    return aus


def pruefe_befund(werte) -> None:
    """Die Aussage des Kapitels muss im Bild sichtbar sein."""
    bester_mse = min(werte, key=lambda z: z[1])
    teuerster = max(werte, key=lambda z: z[2])
    if bester_mse[0] != teuerster[0]:
        raise SystemExit(
            f"Predict-then-Optimize-Bild: Das Verfahren mit dem besten MSE "
            f"({bester_mse[0].replace(chr(10), ' ')}) ist nicht das teuerste "
            f"({teuerster[0].replace(chr(10), ' ')}). Das Kapitel behauptet "
            f"genau das - Bild und Text wuerden sich widersprechen.")


def zeichne(werte) -> None:
    figur, achsen = plt.subplots(1, 2, figsize=(9.4, 4.4))
    namen = [n for n, _, _ in werte]
    y = np.arange(len(namen))[::-1]

    for achse, index, titel, einheit, farbe in (
            (achsen[0], 1, "Prognosefehler (MSE)", "", FARBEN["zweit"]),
            (achsen[1], 2, "Kosten je Tag", " €", FARBEN["warnung"])):
        wert = [w[index] for w in werte]
        bester = min(wert)
        farben = [FARBEN["gut"] if abs(w - bester) < 1e-9 else farbe
                  for w in wert]
        achse.barh(y, wert, height=0.62, color=farben, zorder=3)
        for stelle, w in enumerate(wert):
            achse.text(w * 1.01, y[stelle], f"{w:,.2f}{einheit}"
                       .replace(",", "."), va="center", fontsize=9,
                       color=FARBEN["text"])
        achse.set_yticks(y)
        achse.set_yticklabels(namen if index == 1 else [], fontsize=8.5,
                              color=FARBEN["text"])
        achse.set_xlim(0, max(wert) * 1.22)
        achse.set_title(titel, fontsize=10.5, color=FARBEN["text"])
        achse.grid(axis="y", visible=False)

    figur.suptitle("Grün ist jeweils der beste Wert — und es ist nicht "
                   "dasselbe Verfahren",
                   fontsize=11, color=FARBEN["text"], y=1.02)
    speichere(figur, "kap_prognose_entscheidungsguete")


if __name__ == "__main__":
    werte = rechne()
    for name, mse, kosten in werte:
        print(f"  {name.replace(chr(10), ' '):<40} MSE {mse:>8.2f}   "
              f"Kosten {kosten:>7.2f} EUR")
    pruefe_befund(werte)
    print("  Befund bestätigt: bester MSE = höchste Kosten.")
    zeichne(werte)