Kapitel 23: Testen, Messen, Ausliefern

assert steuer(50_000) == 12_345.60

#!/usr/bin/env python3

# test_or_kern.py
"""
Kapitel Testen: Eine Testsuite fuer ein Optimierungsmodell.

Das Grundproblem beim Testen von Optimierungsmodellen: Man kennt die richtige
Antwort nicht. Bei einer Funktion `steuer(brutto)` schreibt man den erwarteten
Wert hin. Bei "der beste Produktionsplan fuer 400 Auftraege" gibt es niemanden,
der ihn unabhaengig ausrechnen koennte - sonst braeuchte man den Solver nicht.

Deshalb testet man nicht den WERT, sondern vier andere Dinge:

  1. EIGENSCHAFTEN  Die Loesung haelt jede Nebenbedingung ein, und der
                    ausgewiesene Zielwert passt zu den Mengen.
  2. INVARIANTEN    Was das Ergebnis NICHT aendern darf: die Reihenfolge der
                    Produkte, die Waehrungseinheit, ein zusaetzliches Produkt
                    mit Deckungsbeitrag 0.
  3. REGRESSION     Eine kleine Instanz mit von Hand nachgerechnetem Optimum.
  4. FEHLERFAELLE   Unsinnige Eingaben werden abgewiesen, unloesbare Modelle
                    als unloesbar erkannt - und die Abnahmepruefung schlaegt
                    tatsaechlich an, wenn man ihr eine kaputte Loesung gibt.

Der vierte Punkt ist der wichtigste und wird am haeufigsten vergessen: Eine
Pruefung, die noch nie etwas gefunden hat, ist keine Pruefung, sondern eine
Vermutung. Die Tests unten fuettern `pruefe_loesung()` deshalb absichtlich mit
falschen Loesungen und verlangen, dass sie anschlaegt.

Aufruf:
    pytest test_or_kern.py -v
    python3 test_or_kern.py          # ruft pytest selbst auf

Benoetigt: pytest, numpy, pydantic, ortools, scipy (ueber or_kern)
"""

from __future__ import annotations

import pytest

from or_kern import (Loesung, Produkt, Produktionsproblem, SolverStatus,
                     loese_mit_glop, loese_mit_scipy, pruefe_loesung,
                     status_von_scipy)


# --- Die Instanz, gegen die getestet wird ----------------------------------

def schreinerei() -> Produktionsproblem:
    """Der Fall aus Kapitel Einfuehrung - klein genug fuer die Handrechnung.

    Optimum: 30 Tische, 60 Stuehle, Deckungsbeitrag 10.800 EUR.
    Beide Ressourcen sind voll ausgelastet (3*30 + 1*60 = 150,
    6*30 + 1*60 = 240).
    """
    return Produktionsproblem(
        produkte=[
            Produkt(name="Tisch", deckungsbeitrag=240.0,
                    verbrauch={"Montagestunden": 3.0, "Plattenmaterial": 6.0}),
            Produkt(name="Stuhl", deckungsbeitrag=60.0,
                    verbrauch={"Montagestunden": 1.0, "Plattenmaterial": 1.0}),
        ],
        kapazitaeten={"Montagestunden": 150.0, "Plattenmaterial": 240.0},
    )


# Beide Modellbauer durchlaufen JEDEN Test. Ein Test, der nur mit einem Solver
# laeuft, prueft die Bibliothek mit - nicht das Modell.
MODELLBAUER = [loese_mit_glop, loese_mit_scipy]
NAMEN = ["glop", "scipy"]


@pytest.fixture(params=MODELLBAUER, ids=NAMEN)
def bauer(request):
    return request.param


# --- 1. Eigenschaften ------------------------------------------------------

def test_loesung_haelt_alle_nebenbedingungen_ein(bauer):
    """Die Abnahmepruefung darf nichts zu beanstanden haben."""
    problem = schreinerei()
    loesung = bauer(problem)
    assert loesung.status is SolverStatus.OPTIMAL
    assert pruefe_loesung(problem, loesung) == []


def test_zielwert_passt_zu_den_mengen(bauer):
    """Der ausgewiesene Zielwert wird unabhaengig nachgerechnet.

    Klingt trivial, ist es nicht: Wer die Zielfunktion im Modell anders
    zusammensetzt als im Bericht (etwa Gebuehren einmal abgezogen, einmal
    nicht), merkt es sonst nie.
    """
    problem = schreinerei()
    loesung = bauer(problem)
    nachgerechnet = sum(p.deckungsbeitrag * loesung.werte[p.name]
                        for p in problem.produkte)
    assert loesung.zielwert == pytest.approx(nachgerechnet, abs=1e-6)


def test_beide_solver_liefern_dasselbe():
    """Der Kern der Architektur aus Kapitel Praxisfallen, als Test."""
    problem = schreinerei()
    glop, scipy_ = loese_mit_glop(problem), loese_mit_scipy(problem)
    assert glop.zielwert == pytest.approx(scipy_.zielwert, abs=1e-6)
    for ressource in problem.ressourcen:
        assert glop.schattenpreise[ressource] == pytest.approx(
            scipy_.schattenpreise[ressource], abs=1e-6)


# --- 2. Invarianten --------------------------------------------------------

def test_produktreihenfolge_aendert_nichts(bauer):
    """Dieselben Daten in anderer Zeilenreihenfolge - gleiches Ergebnis.

    Das faengt die Spaltenvertauschungsfalle aus Kapitel Finanzdaten ab: ein
    Modell, das Positionen statt Namen benutzt, faellt hier durch.
    """
    problem = schreinerei()
    gedreht = Produktionsproblem(produkte=list(reversed(problem.produkte)),
                                 kapazitaeten=problem.kapazitaeten)
    assert bauer(problem).zielwert == pytest.approx(bauer(gedreht).zielwert)


def test_waehrungseinheit_skaliert_linear(bauer):
    """Deckungsbeitraege in Cent statt Euro: Zielwert mal 100, Mengen gleich.

    Der Test prueft nicht die Mathematik - die ist offensichtlich -, sondern
    die NUMERIK. Wer schlecht skalierte Modelle baut (Kapitel Fundament),
    bekommt hier Abweichungen weit ueber der Toleranz.
    """
    problem = schreinerei()
    in_cent = Produktionsproblem(
        produkte=[Produkt(name=p.name, deckungsbeitrag=p.deckungsbeitrag * 100,
                          verbrauch=p.verbrauch) for p in problem.produkte],
        kapazitaeten=problem.kapazitaeten)
    basis, skaliert = bauer(problem), bauer(in_cent)
    assert skaliert.zielwert == pytest.approx(basis.zielwert * 100, rel=1e-9)
    for p in problem.produkte:
        assert skaliert.werte[p.name] == pytest.approx(basis.werte[p.name],
                                                       abs=1e-6)


def test_produkt_ohne_deckungsbeitrag_aendert_das_optimum_nicht(bauer):
    """Ein Produkt, das nichts einbringt, darf den Zielwert nicht heben."""
    problem = schreinerei()
    mit_nullprodukt = Produktionsproblem(
        produkte=problem.produkte + [
            Produkt(name="Muster", deckungsbeitrag=0.0,
                    verbrauch={"Montagestunden": 1.0, "Plattenmaterial": 1.0})],
        kapazitaeten=problem.kapazitaeten)
    assert bauer(mit_nullprodukt).zielwert <= bauer(problem).zielwert + 1e-6


# --- 3. Regression ---------------------------------------------------------

def test_bekanntes_optimum(bauer):
    """Von Hand nachgerechnet (Kapitel Einfuehrung, Handrechnung).

    Das ist der einzige Test, der einen ZAHLENWERT festschreibt - und er darf
    es, weil dieser Wert von Hand belegt ist. Fuer grosse Instanzen gibt es
    diesen Test nicht; dort tragen die Eigenschaften und Invarianten.
    """
    loesung = bauer(schreinerei())
    assert loesung.zielwert == pytest.approx(10_800.0, abs=1e-6)
    assert loesung.werte["Tisch"] == pytest.approx(30.0, abs=1e-6)
    assert loesung.werte["Stuhl"] == pytest.approx(60.0, abs=1e-6)


def test_schattenpreise_bekannt(bauer):
    """40 EUR je Montagestunde, 20 EUR je Einheit Plattenmaterial."""
    loesung = bauer(schreinerei())
    assert loesung.schattenpreise["Montagestunden"] == pytest.approx(40.0, abs=1e-6)
    assert loesung.schattenpreise["Plattenmaterial"] == pytest.approx(20.0, abs=1e-6)


# --- 4. Fehlerfaelle: schlaegt die Pruefung ueberhaupt an? -----------------

def test_pruefung_findet_kapazitaetsverletzung():
    """Eine Loesung, die zu viel verbraucht, MUSS beanstandet werden."""
    problem = schreinerei()
    kaputt = Loesung(status=SolverStatus.OPTIMAL,
                     werte={"Tisch": 50.0, "Stuhl": 60.0},   # 3*50+60 = 210 > 150
                     zielwert=50 * 240.0 + 60 * 60.0)
    beanstandungen = pruefe_loesung(problem, kaputt)
    assert any("Montagestunden" in b for b in beanstandungen)


def test_pruefung_findet_falschen_zielwert():
    """Mengen und Zielwert passen nicht zusammen - der haeufigste stille Fehler."""
    problem = schreinerei()
    kaputt = Loesung(status=SolverStatus.OPTIMAL,
                     werte={"Tisch": 30.0, "Stuhl": 60.0},
                     zielwert=99_999.0)
    assert any("Zielwert" in b for b in pruefe_loesung(problem, kaputt))


def test_pruefung_findet_negative_mengen():
    problem = schreinerei()
    kaputt = Loesung(status=SolverStatus.OPTIMAL,
                     werte={"Tisch": -5.0, "Stuhl": 60.0},
                     zielwert=-5 * 240.0 + 60 * 60.0)
    assert any("negativ" in b for b in pruefe_loesung(problem, kaputt))


def test_pruefung_findet_gebrochene_ganzzahligkeit():
    """Die Falle aus Kapitel MILP: 0,99999998 ist nicht 1."""
    problem = schreinerei()
    fast_ganz = Loesung(status=SolverStatus.OPTIMAL,
                        werte={"Tisch": 29.4, "Stuhl": 61.8},
                        zielwert=29.4 * 240.0 + 61.8 * 60.0)
    beanstandungen = pruefe_loesung(problem, fast_ganz, ganzzahlig=["Tisch"])
    assert any("ganzzahlig" in b for b in beanstandungen)


def test_pruefung_verweigert_unbrauchbaren_status():
    """Ohne verwertbares Ergebnis wird gar nicht erst gerechnet."""
    problem = schreinerei()
    ohne = Loesung(status=SolverStatus.UNZULAESSIG)
    beanstandungen = pruefe_loesung(problem, ohne)
    assert len(beanstandungen) == 1 and "unzulaessig" in beanstandungen[0]


@pytest.mark.parametrize("daten, stichwort", [
    (dict(produkte=[Produkt(name="Tisch", deckungsbeitrag=240.0,
                            verbrauch={"Montagestunden": 3.0})],
          kapazitaeten={"Montagestunden": 0.0}), "greater than"),
    (dict(produkte=[Produkt(name="Regal", deckungsbeitrag=130.0,
                            verbrauch={"Lackieren": 2.0})],
          kapazitaeten={"Montagestunden": 150.0}), "Kapazitaetsangabe"),
    (dict(produkte=[], kapazitaeten={"Montagestunden": 150.0}), "at least 1"),
])
def test_unsinnige_eingaben_werden_beim_einlesen_abgewiesen(daten, stichwort):
    """Drei Fehler, die NICHT erst beim Loesen auffallen duerfen."""
    with pytest.raises(Exception) as fehler:
        Produktionsproblem(**daten)
    assert stichwort in str(fehler.value)


def test_doppelter_produktname_wird_abgewiesen():
    produkt = Produkt(name="Tisch", deckungsbeitrag=240.0,
                      verbrauch={"Montagestunden": 3.0})
    with pytest.raises(Exception) as fehler:
        Produktionsproblem(produkte=[produkt, produkt],
                           kapazitaeten={"Montagestunden": 150.0})
    assert "mehrfach" in str(fehler.value)


@pytest.mark.parametrize("rohstatus, erwartet", [
    (0, SolverStatus.OPTIMAL),
    (1, SolverStatus.ZEITLIMIT),
    (2, SolverStatus.UNZULAESSIG),
    (3, SolverStatus.UNBESCHRAENKT),
    (4, SolverStatus.FEHLERHAFT),
    (99, SolverStatus.UNBEKANNT),
])
def test_statusuebersetzung_scipy(rohstatus, erwartet):
    """Die Uebersetzertabelle selbst - inklusive des unbekannten Falls."""
    class Ergebnis:
        status = rohstatus
    assert status_von_scipy(Ergebnis()) is erwartet


def test_status_eigenschaften_sind_konsistent():
    """brauchbar und modellfehler duerfen sich nie ueberschneiden."""
    for status in SolverStatus:
        assert not (status.brauchbar and status.modellfehler)


# Die beiden folgenden Tests gab es zuerst nicht. Sie sind entstanden, weil
# Mutationstest.py zwei eingebaute Fehler UEBERLEBEN liess - siehe den
# Abschnitt "Wer testet die Tests?".

def test_nur_optimal_und_zulaessig_gelten_als_brauchbar():
    """ZEITLIMIT ist NICHT brauchbar - der Solver hat nichts gefunden.

    Der Unterschied entscheidet, ob ein Nachtjob einen Plan ausliefert oder
    Alarm schlaegt. Ein aufgeweichtes 'brauchbar' faellt sonst nirgends auf.
    """
    assert SolverStatus.OPTIMAL.brauchbar
    assert SolverStatus.ZULAESSIG.brauchbar
    for status in (SolverStatus.ZEITLIMIT, SolverStatus.UNZULAESSIG,
                   SolverStatus.UNBESCHRAENKT, SolverStatus.FEHLERHAFT,
                   SolverStatus.UNBEKANNT):
        assert not status.brauchbar, f"{status.value} darf nicht brauchbar sein"


def test_ressourcenreihenfolge_haengt_nicht_an_der_eingabe():
    """Zwei Mappings mit gleichem Inhalt, andere Einfuegereihenfolge.

    Die Verbrauchsmatrix muss zeilenweise identisch sein. Ohne diese
    Zusicherung passt die Matrix irgendwann nicht mehr zum
    Kapazitaetsvektor - die Vertauschungsfalle aus Kapitel Finanzdaten,
    nur eine Ebene tiefer.
    """
    produkte = schreinerei().produkte
    vorwaerts = Produktionsproblem(
        produkte=produkte,
        kapazitaeten={"Montagestunden": 150.0, "Plattenmaterial": 240.0})
    rueckwaerts = Produktionsproblem(
        produkte=produkte,
        kapazitaeten={"Plattenmaterial": 240.0, "Montagestunden": 150.0})
    assert vorwaerts.ressourcen == rueckwaerts.ressourcen
    assert (vorwaerts.verbrauchsmatrix() == rueckwaerts.verbrauchsmatrix()).all()
    assert (vorwaerts.kapazitaetsvektor() == rueckwaerts.kapazitaetsvektor()).all()


if __name__ == "__main__":
    import sys
    sys.exit(pytest.main([__file__, "-v", "--tb=short", "-p", "no:cacheprovider"]))

#!/usr/bin/env python3

# Mutationstest.py
"""
Kapitel Testen: Wer testet die Tests?

`test_or_kern.py` besteht aus 31 Tests, und alle sind gruen. Das ist noch kein
Ergebnis. Eine Testsuite, die nie rot war, ist keine Pruefung, sondern eine
Vermutung - sie koennte aus lauter `assert True` bestehen und saehe genauso aus.

Dieses Programm prueft die Pruefung. Es baut nacheinander KLEINE, gezielte
Fehler in `or_kern.py` ein - jeder einzelne ist ein Fehler, den man beim
Programmieren tatsaechlich macht - und laesst die Testsuite darauf los. Fuer
jede Mutation gibt es genau zwei moegliche Ausgaenge:

    GETOETET    Mindestens ein Test wird rot. Die Suite haette diesen Fehler
                in der Praxis gefunden.
    UEBERLEBT   Alle Tests bleiben gruen. Die Suite haette diesen Fehler
                DURCHGEHEN LASSEN - hier fehlt ein Test.

Ueberlebende Mutationen sind das eigentliche Ergebnis: Sie zeigen die Luecken,
und zwar ohne dass man sie erraten muss.

Das Verfahren heisst Mutationstest und ist alt; neu ist nur, wie gut es zu
Optimierungsmodellen passt. Dort kann man den Zielwert meist nicht gegen einen
bekannten Sollwert pruefen - also weiss man ohne so ein Experiment nicht, ob
die Ersatzpruefungen (Eigenschaften, Invarianten) ueberhaupt greifen.

WICHTIG: Es wird nichts veraendert. Jede Mutation laeuft in einem eigenen
temporaeren Verzeichnis mit einer Kopie von or_kern.py.

Benoetigt: pytest; im selben Verzeichnis or_kern.py und test_or_kern.py
"""

from __future__ import annotations

import os
import shutil
import subprocess
import sys
import tempfile

# (Name, gesuchter Text, Ersatz) - jede Zeile ein realistischer Fehler.
MUTATIONEN = [
    ("Nebenbedingung umgedreht (<= wird >=)",
     "<= problem.kapazitaeten[ressource], name=ressource)",
     ">= problem.kapazitaeten[ressource], name=ressource)"),

    ("Maximierung wird Minimierung",
     "solver.Maximize(sum(menge[p.name] * p.deckungsbeitrag",
     "solver.Minimize(sum(menge[p.name] * p.deckungsbeitrag"),

    ("Vorzeichen der Dualwerte vergessen",
     "schattenpreise={r: float(-m) for r, m in",
     "schattenpreise={r: float(m) for r, m in"),

    ("Abnahmepruefung aufgeweicht (Toleranz 1e-6 -> 1.0)",
     "toleranz: float = 1e-6) -> list[str]:",
     "toleranz: float = 1.0) -> list[str]:"),

    ("Kapazitaetspruefung uebersprungen",
     "if ist > grenze + toleranz:",
     "if False:"),

    ("Pruefung auf doppelte Produktnamen entfernt",
     "if len(set(namen)) != len(namen):",
     "if False:"),

    ("Ganzzahligkeitspruefung entfernt",
     "if abs(wert - round(wert)) > toleranz:",
     "if False:"),

    ("Zielwertabgleich entfernt",
     "if abs(nachgerechnet - loesung.zielwert) > toleranz * max(1.0, abs(nachgerechnet)):",
     "if False:"),

    ("ZEITLIMIT gilt faelschlich als brauchbar",
     "return self in (SolverStatus.OPTIMAL, SolverStatus.ZULAESSIG)",
     "return self in (SolverStatus.OPTIMAL, SolverStatus.ZULAESSIG,\n"
     "                        SolverStatus.ZEITLIMIT)"),

    ("Ressourcenreihenfolge nicht mehr stabil",
     "return sorted(self.kapazitaeten)",
     "return list(self.kapazitaeten)"),
]


def fuehre_suite_aus(verzeichnis: str) -> tuple[int, int, str]:
    """Laesst pytest im Verzeichnis laufen; liefert (bestanden, fehlgeschlagen)."""
    ergebnis = subprocess.run(
        [sys.executable, "-m", "pytest", "test_or_kern.py", "-q",
         "--tb=no", "-p", "no:cacheprovider"],
        cwd=verzeichnis, capture_output=True, text=True, timeout=600)
    letzte = [z for z in ergebnis.stdout.strip().splitlines() if z.strip()]
    zeile = letzte[-1] if letzte else ""
    bestanden = fehlgeschlagen = 0
    for teil in zeile.replace("=", " ").split(","):
        for wort in teil.split():
            if wort.isdigit():
                zahl = int(wort)
            elif wort.startswith("passed"):
                bestanden = zahl
            elif wort.startswith("failed") or wort.startswith("error"):
                fehlgeschlagen += zahl
    return bestanden, fehlgeschlagen, zeile


def mutiere_und_pruefe(quelle: str, ziel: str, suche: str,
                       ersatz: str) -> tuple[int, int, str] | None:
    """Legt eine mutierte Kopie an und laesst die Suite darauf laufen."""
    with tempfile.TemporaryDirectory() as verzeichnis:
        text = open(quelle, encoding="utf-8").read()
        if text.count(suche) != 1:
            return None                      # Muster passt nicht (mehr)
        open(os.path.join(verzeichnis, "or_kern.py"), "w",
             encoding="utf-8").write(text.replace(suche, ersatz))
        shutil.copy(ziel, verzeichnis)
        return fuehre_suite_aus(verzeichnis)


if __name__ == "__main__":
    hier = os.path.dirname(os.path.abspath(__file__))
    quelle = os.path.join(hier, "or_kern.py")
    suite = os.path.join(hier, "test_or_kern.py")
    for pfad in (quelle, suite):
        if not os.path.exists(pfad):
            raise SystemExit(f"Nicht gefunden: {pfad}\n"
                             f"Beide Dateien muessen im selben Verzeichnis liegen.")

    print("=" * 84)
    print("  MUTATIONSTEST: WUERDE DIE SUITE DIESE FEHLER FINDEN?")
    print("=" * 84)

    with tempfile.TemporaryDirectory() as v:
        shutil.copy(quelle, v)
        shutil.copy(suite, v)
        bestanden, fehlgeschlagen, _ = fuehre_suite_aus(v)
    print(f"Ausgangslage: {bestanden} Tests, {fehlgeschlagen} rot.\n")
    if fehlgeschlagen:
        raise SystemExit("Die Suite ist schon ohne Mutation rot - erst das reparieren.")

    print(f"{'eingebauter Fehler':<48} {'rote Tests':>11} {'Urteil':>12}")
    print("-" * 84)

    getoetet, ueberlebt = [], []
    for name, suche, ersatz in MUTATIONEN:
        ergebnis = mutiere_und_pruefe(quelle, suite, suche, ersatz)
        if ergebnis is None:
            print(f"{name:<48} {'-':>11} {'nicht anwendbar':>12}")
            continue
        _, rot, _ = ergebnis
        urteil = "GETOETET" if rot else "UEBERLEBT"
        (getoetet if rot else ueberlebt).append(name)
        print(f"{name:<48} {rot:>11} {urteil:>12}")

    print("-" * 84)
    quote = len(getoetet) / max(1, len(getoetet) + len(ueberlebt)) * 100
    print(f"Getoetet: {len(getoetet)} von {len(getoetet) + len(ueberlebt)} "
          f"({quote:.0f} %)")

    print("\n" + "=" * 84)
    print("  WAS DAS ERGEBNIS BEDEUTET")
    print("=" * 84)
    if ueberlebt:
        print("Diese Fehler haette die Suite DURCHGEHEN LASSEN:\n")
        for name in ueberlebt:
            print(f"  * {name}")
        print("\nJeder ueberlebende Eintrag ist eine Testluecke - und zwar eine")
        print("gefundene, keine vermutete. Das ist der ganze Zweck des Verfahrens:")
        print("Es sagt einem, welchen Test man als naechstes schreiben sollte,")
        print("statt dass man raten muss.")
    else:
        print("Alle eingebauten Fehler wurden gefunden. Das heisst nicht, dass die")
        print("Suite vollstaendig ist - nur, dass sie diese zehn Fehler faengt.")
    print()
    print("Die Quote selbst ist keine Kennzahl fuer ein Dashboard. Zehn von Hand")
    print("gewaehlte Mutationen sind keine Stichprobe aus der Menge aller")
    print("moeglichen Fehler. Was zaehlt, ist die LISTE der Ueberlebenden.")
    print("=" * 84)

#!/usr/bin/env python3

# Benchmark_Skalierung.py
"""
Kapitel Testen: Ein Vergleich, dem man glauben kann.

Solververgleiche stehen in jedem Blog, und die meisten sind wertlos - nicht
weil falsch gemessen wurde, sondern weil zu wenig dazugesagt wird. Dieses
Programm misst dasselbe Transportproblem in drei Groessen mit vier
Bibliotheken und haelt sich dabei an fuenf Regeln, die den Unterschied machen:

  1. EIGENER PROZESS je Bibliothek. Nicht nur wegen des Importkonflikts
     zwischen ortools und highspy (Kapitel Oekosystem) - auch, damit der
     Speicherverbrauch der einen nicht in der Messung der anderen auftaucht.
  2. AUFBAU UND LOESEN GETRENNT messen. Bei grossen Instanzen ist der Aufbau
     des Modells in Python regelmaessig teurer als das Loesen. Wer nur die
     Gesamtzeit misst, optimiert am Ende die falsche Haelfte.
  3. ZIELWERTE GEGENEINANDER PRUEFEN. Eine Bibliothek, die schneller ist und
     etwas anderes ausrechnet, hat den Vergleich nicht gewonnen. Diese Pruefung
     ist der wichtigste Teil des Programms.
  4. SPEICHER MITMESSEN. Bei 10.000 Variablen entscheidet oft er und nicht die
     Zeit darueber, was auf einer Maschine noch laeuft.
  5. DIESELBE INSTANZ fuer alle. Feste Saat, kein Neuwuerfeln zwischendurch.

Und die Einschraenkung, die dazugehoert: Gemessen wird EIN Problemtyp in EINER
Formulierung auf EINER Maschine. Das Ergebnis ist keine Rangliste der Solver,
sondern eine Entscheidungshilfe fuer genau diesen Fall. Wer es verallgemeinert,
macht denselben Fehler wie jemand, der aus einem Backtest auf die Zukunft
schliesst (Kapitel Handelsmaschine).

Benoetigt: numpy; in den Kindprozessen scipy, highspy, ortools, cvxpy
"""

from __future__ import annotations

import multiprocessing
import resource
import time
from concurrent.futures import ProcessPoolExecutor

import numpy as np

GROESSEN = [(10, 10), (32, 32), (100, 100)]   # (Lager, Kunden) -> 100 / 1.024 / 10.000 Variablen


# Instanz und Speichermessung stehen als gewoehnliche Funktionen hier - nicht
# in einem String, den ein Kindprozess ausfuehrt. Jede Messfunktion baut die
# Instanz aus derselben Saat neu auf, damit ueber die Prozessgrenze nichts
# reist, was das Ergebnis verfaelschen koennte.

def instanz(m: int, n: int):
    rng = np.random.default_rng(20)
    kosten = rng.integers(5, 95, (m, n)).astype(float)
    angebot = rng.integers(50, 150, m).astype(float)
    bedarf = angebot.sum() * rng.dirichlet(np.ones(n))
    return kosten, angebot, bedarf


def speicher_mb() -> float:
    # ru_maxrss ist unter Linux in Kilobyte. Gemessen wird der Kindprozess -
    # deshalb muss jede Messung einen eigenen bekommen.
    return resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / 1024


def messe_scipy(m: int, n: int):
    from scipy.optimize import linprog
    kosten, angebot, bedarf = instanz(m, n)
    t0 = time.perf_counter()
    c = kosten.reshape(-1)
    A_ub = np.zeros((m, m * n)); A_eq = np.zeros((n, m * n))
    for i in range(m):
        A_ub[i, i * n:(i + 1) * n] = 1.0
    for j in range(n):
        A_eq[j, j::n] = 1.0
    aufbau = time.perf_counter() - t0
    t0 = time.perf_counter()
    r = linprog(c=c, A_ub=A_ub, b_ub=angebot, A_eq=A_eq, b_eq=bedarf,
                bounds=(0, None), method="highs")
    loesen = time.perf_counter() - t0
    return float(r.fun), aufbau, loesen, speicher_mb()


def messe_highspy(m: int, n: int):
    import highspy
    kosten, angebot, bedarf = instanz(m, n)
    t0 = time.perf_counter()
    h = highspy.Highs(); h.setOptionValue("output_flag", False)
    h.addVars(m * n, np.zeros(m * n), np.full(m * n, highspy.kHighsInf))
    for k in range(m * n):
        h.changeColCost(k, float(kosten.reshape(-1)[k]))
    for i in range(m):
        idx = np.arange(i * n, (i + 1) * n, dtype=np.int32)
        h.addRow(-highspy.kHighsInf, float(angebot[i]), n, idx, np.ones(n))
    for j in range(n):
        idx = np.arange(j, m * n, n, dtype=np.int32)
        h.addRow(float(bedarf[j]), float(bedarf[j]), m, idx, np.ones(m))
    aufbau = time.perf_counter() - t0
    t0 = time.perf_counter(); h.run(); loesen = time.perf_counter() - t0
    return h.getInfo().objective_function_value, aufbau, loesen, speicher_mb()


def messe_ortools(m: int, n: int):
    from ortools.linear_solver import pywraplp
    kosten, angebot, bedarf = instanz(m, n)
    t0 = time.perf_counter()
    s = pywraplp.Solver.CreateSolver("GLOP")
    x = [[s.NumVar(0, s.infinity(), f"x{i}_{j}") for j in range(n)]
         for i in range(m)]
    for i in range(m):
        s.Add(sum(x[i]) <= float(angebot[i]))
    for j in range(n):
        s.Add(sum(x[i][j] for i in range(m)) == float(bedarf[j]))
    s.Minimize(sum(float(kosten[i, j]) * x[i][j]
                   for i in range(m) for j in range(n)))
    aufbau = time.perf_counter() - t0
    t0 = time.perf_counter(); s.Solve(); loesen = time.perf_counter() - t0
    return s.Objective().Value(), aufbau, loesen, speicher_mb()


def messe_cvxpy(m: int, n: int):
    import cvxpy as cp
    kosten, angebot, bedarf = instanz(m, n)
    t0 = time.perf_counter()
    x = cp.Variable((m, n), nonneg=True)
    problem = cp.Problem(cp.Minimize(cp.sum(cp.multiply(kosten, x))),
                         [cp.sum(x, axis=1) <= angebot,
                          cp.sum(x, axis=0) == bedarf])
    aufbau = time.perf_counter() - t0
    t0 = time.perf_counter(); problem.solve(); loesen = time.perf_counter() - t0
    return float(problem.value), aufbau, loesen, speicher_mb()


ANSAETZE = {"scipy.linprog": messe_scipy, "highspy": messe_highspy,
            "ortools/GLOP": messe_ortools, "cvxpy": messe_cvxpy}


def messe(funktion, m: int, n: int):
    """Fuehrt eine Messfunktion in einem FRISCHEN Prozess aus.

    'spawn' und max_tasks_per_child=1 zusammen garantieren, was Regel 1
    verlangt: Jede Messung sieht einen leeren Interpreter. Ohne das
    zweite wuerde der Pool seinen Arbeiter wiederverwenden - dann waere
    der Speicherwert der zweiten Bibliothek um die erste zu hoch, und
    ortools und highspy saessen im selben Prozess.
    """
    with ProcessPoolExecutor(
            max_workers=1,
            mp_context=multiprocessing.get_context("spawn"),
            max_tasks_per_child=1) as pool:
        try:
            return pool.submit(funktion, m, n).result(timeout=600), None
        except Exception as fehler:
            return None, str(fehler).strip().splitlines()[-1][:60]


if __name__ == "__main__":
    print("=" * 92)
    print("  SKALIERUNGSVERGLEICH: TRANSPORTPROBLEM, VIER BIBLIOTHEKEN")
    print("=" * 92)
    print("Jede Zeile ein eigener Prozess. Zeiten und Speicher sind "
          "hardwareabhaengig,")
    print("die Zielwerte und ihr Verhaeltnis zueinander nicht.\n")

    for m, n in GROESSEN:
        kopf = f"--- {m} Lager x {n} Kunden = {m * n:,} Variablen "
        print(kopf + "-" * max(3, 92 - len(kopf)))
        print(f"  {'Bibliothek':<16} {'Zielwert':>14} {'Aufbau':>9} "
              f"{'Loesen':>9} {'Anteil':>8} {'Speicher':>10}")
        print("  " + "-" * 72)
        zielwerte = {}
        for name, funktion in ANSAETZE.items():
            werte, fehler = messe(funktion, m, n)
            if werte is None:
                print(f"  {name:<16} nicht verfuegbar: {fehler}")
                continue
            ziel, aufbau, loesen, speicher = werte
            zielwerte[name] = ziel
            anteil = aufbau / (aufbau + loesen) * 100
            print(f"  {name:<16} {ziel:>14,.2f} {aufbau:>8.3f}s "
                  f"{loesen:>8.3f}s {anteil:>7.0f}% {speicher:>9.0f} MB")

        # Die wichtigste Zeile: Rechnen alle dasselbe aus?
        spanne = max(zielwerte.values()) - min(zielwerte.values())
        bezug = max(abs(v) for v in zielwerte.values())
        print(f"  {'':16} Spannweite der Zielwerte: {spanne:.2e} "
              f"(relativ {spanne / bezug:.1e})")
        if spanne / bezug > 1e-6:
            print("  ACHTUNG: Die Bibliotheken widersprechen sich - "
                  "der Zeitvergleich ist wertlos.")
        print()

    print("=" * 92)
    print("  WAS MAN AUS SO EINER TABELLE ABLESEN DARF - UND WAS NICHT")
    print("=" * 92)
    print("DARF man ablesen:")
    print("  * Die Spalte 'Anteil' - wie viel der Zeit in den AUFBAU geht statt")
    print("    ins Loesen. Wenn dort 80 % stehen, ist ein schnellerer Solver die")
    print("    falsche Antwort; dann gehoert das Modell vektorisiert aufgebaut")
    print("    (Kapitel Oekosystem).")
    print("  * Die Groessenordnung des Speicherbedarfs. Sie entscheidet, was auf")
    print("    einer bestimmten Maschine ueberhaupt laeuft.")
    print("  * Wie sich beides mit der Groesse ENTWICKELT. Der Trend ist")
    print("    uebertragbarer als der Absolutwert.")
    print()
    print("NICHT ablesen darf man:")
    print("  * 'Bibliothek X ist schneller als Y.' Gemessen wurde EIN")
    print("    Problemtyp in EINER Formulierung. Ein MILP, ein QP oder eine")
    print("    andere Modellierung desselben Problems koennen die Reihenfolge")
    print("    umdrehen.")
    print("  * Etwas ueber Ihre Maschine. Diese Zahlen stammen von einer")
    print("    anderen. Der Sinn des Programms ist, dass Sie es auf Ihrer")
    print("    laufen lassen.")
    print("=" * 92)

#!/usr/bin/env python3

# Optimierungsdienst.py
"""
Kapitel Testen: Das Modell als Dienst - und warum es kein normaler Endpunkt ist.

Ein Optimierungsmodell hinter eine HTTP-Schnittstelle zu haengen sieht nach
einer Fingeruebung aus. Es gibt aber einen Unterschied, der alles bestimmt:

    Eine gewoehnliche Anfrage dauert Millisekunden.
    Eine Optimierung dauert Sekunden bis Minuten - und manchmal laenger,
    als jemand warten will.

Damit scheidet die naheliegende Bauform aus. Wer den Solver direkt im
Request-Handler aufruft, baut einen Dienst, der bei der ersten grossen Instanz
in einen Timeout des Reverse Proxy laeuft, und bei zehn gleichzeitigen Anfragen
alle Arbeiter blockiert.

Die tragfaehige Bauform ist deshalb ZWEISTUFIG:

    POST /plaene         nimmt den Auftrag an, prueft ihn, gibt sofort
                         202 Accepted und eine Auftragsnummer zurueck
    GET  /plaene/{id}    sagt, wie weit es ist - und liefert am Ende die
                         Loesung samt Status und Gap

Drei Dinge, die dabei leicht untergehen und hier bewusst sichtbar sind:

  1. Das Domaenenmodell aus or_kern.py ist zugleich das API-Schema. Eine
     unsinnige Anfrage wird von FastAPI mit 422 abgewiesen, bevor irgendein
     Solver startet - dieselbe Pruefung wie beim Excel-Import.
  2. Die Antwort enthaelt STATUS und GAP, nicht nur Zahlen. Ein Aufrufer, der
     nur die Mengen bekommt, kann nicht unterscheiden, ob er ein bewiesenes
     Optimum oder einen Zeitlimit-Abbruch in der Hand haelt.
  3. Jeder Auftrag hat ein ZEITLIMIT. Ohne das belegt eine einzige unguenstige
     Instanz einen Arbeiter auf unbestimmte Zeit.

Dieses Programm laeuft ohne Server: Der Selbsttest unten benutzt den
TestClient von FastAPI und spricht die Anwendung direkt an. Fuer den echten
Betrieb steht am Ende, was sich aendert.

Aufruf:
    python3 Optimierungsdienst.py          # Selbsttest, kein Server noetig
    uvicorn Optimierungsdienst:app         # echter Server auf Port 8000

Benoetigt: fastapi, httpx (fuer den Selbsttest), pydantic, ortools (ueber or_kern)
"""

from __future__ import annotations

import time
import uuid
from concurrent.futures import ThreadPoolExecutor
from enum import Enum

from fastapi import FastAPI, HTTPException
from pydantic import BaseModel, Field

from or_kern import (Loesung, Produktionsproblem, loese_mit_glop,
                     pruefe_loesung)

ZEITLIMIT_SEKUNDEN = 30.0
ARBEITER = 2


class Auftragsstand(str, Enum):
    WARTET = "wartet"
    LAEUFT = "laeuft"
    FERTIG = "fertig"
    GESCHEITERT = "gescheitert"


class Auftragsantwort(BaseModel):
    """Was der Aufrufer beim Abholen bekommt.

    Bewusst NICHT nur die Mengen: 'stand' und die Felder aus 'loesung'
    (Status, Gap) sind der Unterschied zwischen einer Zahl und einer
    belastbaren Auskunft.
    """
    id: str
    stand: Auftragsstand
    eingegangen: float
    laufzeit: float | None = None
    loesung: Loesung | None = None
    beanstandungen: list[str] = Field(default_factory=list)
    fehler: str | None = None


app = FastAPI(title="Optimierungsdienst",
              summary="Produktionsplanung als Auftrag, nicht als Abfrage")

# Fuer das Buchbeispiel: Auftragsbuch im Speicher, Arbeiter im selben Prozess.
# Was daran im echten Betrieb nicht reicht, steht unten.
AUFTRAEGE: dict[str, Auftragsantwort] = {}
POOL = ThreadPoolExecutor(max_workers=ARBEITER)


def rechne(auftrags_id: str, problem: Produktionsproblem) -> None:
    """Laeuft im Arbeiterthread - nie im Request-Handler."""
    antwort = AUFTRAEGE[auftrags_id]
    antwort.stand = Auftragsstand.LAEUFT
    start = time.perf_counter()
    try:
        loesung = loese_mit_glop(problem)
        antwort.loesung = loesung
        # Dieselbe Abnahmepruefung wie ueberall sonst. Ein Dienst, der sie
        # weglaesst, liefert Fehler schneller aus als ein Mensch sie faende.
        antwort.beanstandungen = pruefe_loesung(problem, loesung)
        antwort.stand = (Auftragsstand.FERTIG if loesung.status.brauchbar
                         and not antwort.beanstandungen
                         else Auftragsstand.GESCHEITERT)
        if not loesung.status.brauchbar:
            antwort.fehler = f"Solverstatus: {loesung.status.value}"
        elif antwort.beanstandungen:
            antwort.fehler = "Abnahmepruefung fehlgeschlagen"
    except Exception as fehler:                      # noqa: BLE001
        antwort.stand = Auftragsstand.GESCHEITERT
        antwort.fehler = f"{type(fehler).__name__}: {fehler}"
    finally:
        antwort.laufzeit = time.perf_counter() - start


@app.post("/plaene", status_code=202, response_model=Auftragsantwort)
def auftrag_annehmen(problem: Produktionsproblem) -> Auftragsantwort:
    """Nimmt an, prueft, gibt sofort zurueck.

    Der Typ 'Produktionsproblem' im Parameter ist der ganze Trick: FastAPI
    validiert die Anfrage damit gegen das Domaenenmodell und antwortet bei
    Unsinn mit 422, ohne dass hier eine Zeile Pruefcode steht.
    """
    auftrags_id = str(uuid.uuid4())
    AUFTRAEGE[auftrags_id] = Auftragsantwort(
        id=auftrags_id, stand=Auftragsstand.WARTET, eingegangen=time.time())
    POOL.submit(rechne, auftrags_id, problem)
    return AUFTRAEGE[auftrags_id]


@app.get("/plaene/{auftrags_id}", response_model=Auftragsantwort)
def auftrag_abholen(auftrags_id: str) -> Auftragsantwort:
    if auftrags_id not in AUFTRAEGE:
        raise HTTPException(status_code=404, detail="Unbekannter Auftrag")
    return AUFTRAEGE[auftrags_id]


@app.get("/gesundheit")
def gesundheit() -> dict[str, object]:
    """Was ein Ueberwachungssystem abfragt - siehe Betriebsueberwachung.py."""
    offen = sum(1 for a in AUFTRAEGE.values()
                if a.stand in (Auftragsstand.WARTET, Auftragsstand.LAEUFT))
    return {"zustand": "bereit", "auftraege_gesamt": len(AUFTRAEGE),
            "offen": offen, "arbeiter": ARBEITER}


# --- Selbsttest ohne Server -------------------------------------------------

SCHREINEREI = {
    "produkte": [
        {"name": "Tisch", "deckungsbeitrag": 240.0,
         "verbrauch": {"Montagestunden": 3.0, "Plattenmaterial": 6.0}},
        {"name": "Stuhl", "deckungsbeitrag": 60.0,
         "verbrauch": {"Montagestunden": 1.0, "Plattenmaterial": 1.0}},
    ],
    "kapazitaeten": {"Montagestunden": 150.0, "Plattenmaterial": 240.0},
}

DOCKERFILE = """\
FROM python:3.12-slim
WORKDIR /app
COPY requirements.txt .
RUN pip install --no-cache-dir -r requirements.txt
COPY or_kern.py Optimierungsdienst.py .
# Ein Arbeiter je CPU-Kern - der Solver rechnet selbst schon parallel,
# mehr Prozesse machen ihn nicht schneller, sondern langsamer.
CMD ["uvicorn", "Optimierungsdienst:app", "--host", "0.0.0.0", "--port", "8000"]
"""


def warte_auf_ergebnis(klient, auftrags_id: str, grenze: float = 30.0) -> dict:
    ende = time.time() + grenze
    while time.time() < ende:
        antwort = klient.get(f"/plaene/{auftrags_id}").json()
        if antwort["stand"] in ("fertig", "gescheitert"):
            return antwort
        time.sleep(0.02)
    raise TimeoutError("Auftrag wurde nicht fertig")


if __name__ == "__main__":
    from fastapi.testclient import TestClient

    print("=" * 80)
    print("  DAS MODELL ALS DIENST")
    print("=" * 80)

    with TestClient(app) as klient:
        # 1. Gueltiger Auftrag
        print("\n1. Auftrag einreichen")
        angenommen = klient.post("/plaene", json=SCHREINEREI)
        auftrags_id = angenommen.json()["id"]
        print(f"   POST /plaene            -> {angenommen.status_code} "
              f"{'Accepted' if angenommen.status_code == 202 else ''}")
        print(f"   Auftragsnummer          -> UUID, {len(auftrags_id)} Zeichen")
        # Der entscheidende Punkt, und er ist pruefbar: Die Antwort ist da,
        # BEVOR es eine Loesung gibt. Genau dafuer ist 202 gedacht.
        print(f"   Loesung schon dabei?    -> "
              f"{'ja' if angenommen.json()['loesung'] else 'nein - genau so soll es sein'}")

        fertig = warte_auf_ergebnis(klient, auftrags_id)
        loesung = fertig["loesung"]
        print("\n2. Ergebnis abholen")
        print(f"   Stand                   -> {fertig['stand']}")
        print(f"   Solverstatus            -> {loesung['status']}")
        print(f"   Zielwert                -> {loesung['zielwert']:,.2f} EUR")
        print(f"   Mengen                  -> " + ", ".join(
            f"{name} {wert:.0f}" for name, wert in loesung["werte"].items()))
        print(f"   Schattenpreise          -> " + ", ".join(
            f"{name} {wert:.0f}" for name, wert in loesung["schattenpreise"].items()))
        print(f"   Abnahmepruefung         -> "
              f"{'bestanden' if not fertig['beanstandungen'] else fertig['beanstandungen']}")

        # 2. Unsinnige Anfrage - muss abgewiesen werden, BEVOR gerechnet wird
        print("\n3. Unsinnige Anfragen (der Solver startet gar nicht erst)")
        for beschreibung, aenderung in [
                ("Kapazitaet 0", {"kapazitaeten": {"Montagestunden": 0.0}}),
                ("keine Produkte", {"produkte": []}),
                ("Ressource ohne Kapazitaet",
                 {"produkte": [{"name": "Regal", "deckungsbeitrag": 130.0,
                                "verbrauch": {"Lackieren": 2.0}}]})]:
            anfrage = {**SCHREINEREI, **aenderung}
            antwort = klient.post("/plaene", json=anfrage)
            print(f"   {beschreibung:<24} -> {antwort.status_code} "
                  f"{'Unprocessable Content' if antwort.status_code == 422 else ''}")

        print("\n4. Unbekannter Auftrag")
        print(f"   GET /plaene/gibtsnicht  -> "
              f"{klient.get('/plaene/gibtsnicht').status_code} Not Found")

        print("\n5. Gesundheitsabfrage")
        print(f"   GET /gesundheit         -> {klient.get('/gesundheit').json()}")

    # --- Warum ein Thread je Auftrag genuegt -----------------------------
    print("\n" + "-" * 80)
    print("6. Warum Threads hier reichen (und wann nicht)\n")
    from concurrent.futures import ThreadPoolExecutor as Pool

    # Wichtig: eine Rechnung mit ECHTER Last. Ein Modell, das in zehn
    # Millisekunden fertig ist, misst nur den Aufwand fuers Threadstarten -
    # der erste Entwurf dieser Messung ist genau daran gescheitert und zeigte
    # eine Verlangsamung, wo in Wirklichkeit eine Beschleunigung steht.
    def eine_rechnung() -> None:
        import numpy as np
        from ortools.linear_solver import pywraplp
        rng = np.random.default_rng(3)
        n = 260
        solver = pywraplp.Solver.CreateSolver("SCIP")
        x = [solver.IntVar(0, 1, f"x{i}") for i in range(n)]
        gewicht = rng.integers(10, 60, n)
        wert = rng.integers(10, 60, n)
        solver.Add(sum(int(gewicht[i]) * x[i] for i in range(n))
                   <= int(gewicht.sum() * 0.5))
        for _ in range(30):
            auswahl = rng.choice(n, 40, replace=False)
            solver.Add(sum(x[int(i)] for i in auswahl) <= 12)
        solver.Maximize(sum(int(wert[i]) * x[i] for i in range(n)))
        solver.Solve()

    t0 = time.perf_counter()
    for _ in range(4):
        eine_rechnung()
    seriell = time.perf_counter() - t0
    t0 = time.perf_counter()
    with Pool(max_workers=4) as p:
        list(p.map(lambda _: eine_rechnung(), range(4)))
    parallel = time.perf_counter() - t0

    print(f"   vier Rechnungen nacheinander -> {seriell:6.2f} s")
    print(f"   vier Rechnungen in Threads   -> {parallel:6.2f} s "
          f"(Faktor {seriell / parallel:.1f})")
    print()
    if parallel < seriell * 0.75:
        print("   Die Laeufe ueberlappen sich. ortools rechnet in C++ und gibt")
        print("   den GIL waehrend Solve() frei - deshalb genuegt hier ein")
        print("   Threadpool, es braucht keine eigenen Prozesse.")
    else:
        print("   Keine Ueberlappung - hier waeren Prozesse noetig.")
    print()
    print("   ACHTUNG, das gilt nicht allgemein: Eine in reinem Python")
    print("   geschriebene Heuristik (Kapitel Metaheuristiken) haelt den GIL")
    print("   die ganze Zeit. Fuer sie braucht derselbe Dienst einen")
    print("   ProcessPoolExecutor statt eines Threadpools.")

    print("\n" + "=" * 80)
    print("  WAS SICH IM ECHTEN BETRIEB AENDERT")
    print("=" * 80)
    print("Dieses Beispiel haelt das Auftragsbuch im Speicher und rechnet in")
    print("Threads desselben Prozesses. Das reicht zum Zeigen und fuer einen")
    print("einzelnen Rechner - nicht darueber hinaus:")
    print()
    print("  * NEUSTART LOESCHT ALLES. Auftragsbuch in eine Datenbank oder eine")
    print("    Warteschlange (Redis, RabbitMQ), nicht in ein dict.")
    print("  * ZWEI INSTANZEN KENNEN EINANDER NICHT. Sobald der Dienst mehr als")
    print("    einmal laeuft, muss die Warteschlange ausserhalb liegen -")
    print("    typischerweise Celery mit Redis als Vermittler.")
    print("  * DIE WAHL THREADS/PROZESSE HAENGT AM SOLVER - siehe die Messung")
    print("    unter Punkt 6. Sie gehoert gemessen, nicht angenommen.")
    print("  * OHNE ZEITLIMIT KEIN DIENST. Jede Instanz bekommt eines, und der")
    print("    Aufrufer erfaehrt im Status, ob es gegriffen hat.")
    print()
    print("Das zugehoerige Dockerfile ist kurz genug, um es ganz zu zeigen:")
    print()
    for zeile in DOCKERFILE.splitlines():
        print(f"    {zeile}")
    print()
    print("(Es wird hier nicht gebaut - das Buch setzt keine laufende")
    print(" Docker-Installation voraus.)")
    print("=" * 80)