Kapitel 3: Das Python-Ökosystem für OR — Solver, Bindings und Modellierungsschichten

#!/usr/bin/env python3

# Solver_Wahl.py
"""
Kapitel Oekosystem: Entscheidungshilfe zur Solverwahl.
Beantwortet drei Fragen und empfiehlt Bibliothek + Backend.
"""

from dataclasses import dataclass


@dataclass
class Problem:
    diskrete_entscheidungen: bool   # Ja/Nein-Variablen, Zuordnungen, Reihenfolgen?
    zielfunktion: str               # "linear" | "quadratisch" | "konvex" | "beliebig"
    routing: bool = False           # Fahrzeuge, Depots, Zeitfenster?
    groesse: str = "klein"          # "klein" (<10^3 Var.) | "gross"
    wiederholte_laeufe: bool = False


def empfehle(p: Problem) -> tuple[str, str, str]:
    """Gibt (Bibliothek, Backend, Begründung) zurück."""
    if p.routing:
        return ("ortools.constraint_solver (Routing)", "Routing Engine",
                "Spezialisierte Heuristiken für VRP/TSP schlagen jedes selbstgebaute Modell.")

    if p.diskrete_entscheidungen:
        if p.zielfunktion in ("linear",) and not p.routing:
            if p.groesse == "gross":
                return ("highspy", "HiGHS Branch-and-Cut",
                        "MILP mit ökonomischer Struktur (Fixkosten, Kardinalität).")
            return ("ortools.sat (CP-SAT)", "CP-SAT",
                    "Logische Regeln und Zuweisungen: CP-SAT propagiert sehr effizient.")
        return ("ortools.sat (CP-SAT)", "CP-SAT",
                "Diskrete Struktur dominiert; CP-SAT verarbeitet auch nichtlineare Logik.")

    if p.zielfunktion == "linear":
        if p.wiederholte_laeufe or p.groesse == "gross":
            return ("highspy", "HiGHS Dual Simplex",
                    "Modell einmal aufbauen, Parameter ändern, wiederholt lösen.")
        return ("scipy.optimize.linprog", "HiGHS",
                "Kleinstes Setup, keine zusätzliche Abhängigkeit.")

    if p.zielfunktion in ("quadratisch", "konvex"):
        return ("cvxpy", "Clarabel / OSQP",
                "Konvexität wird automatisch geprüft; Portfolio-Standard.")

    return ("scipy.optimize.minimize", "SLSQP / trust-constr",
            "Nicht konvex: nur lokales Optimum, Startpunkt variieren und vergleichen!")


BEISPIELE = {
    "Vertretungsplan Schule":        Problem(True,  "linear", groesse="klein"),
    "Produktionsplanung (LP)":       Problem(False, "linear", groesse="klein"),
    "Produktionsplanung, 50k Var.":  Problem(False, "linear", groesse="gross"),
    "Portfolio Markowitz":           Problem(False, "quadratisch"),
    "Portfolio mit max. 5 Titeln":   Problem(True,  "quadratisch"),
    "Liefertouren mit Zeitfenstern": Problem(True,  "linear", routing=True),
    "Entropie-Allokation (NLP)":     Problem(False, "beliebig"),
    "Backtest, 60x neu optimieren":  Problem(False, "konvex", wiederholte_laeufe=True),
}

if __name__ == "__main__":
    print("=" * 92)
    print("  SOLVER-EMPFEHLUNG")
    print("=" * 92)
    for name, p in BEISPIELE.items():
        lib, backend, grund = empfehle(p)
        print(f"\n{name}")
        print(f"   -> Bibliothek: {lib}")
        print(f"      Backend:    {backend}")
        print(f"      Grund:      {grund}")
    print("\n" + "=" * 92)

#!/usr/bin/env python3

# Ein_System_Vier_Ansaetze.py
"""
Kapitel Oekosystem: Dasselbe LP in vier Bibliotheken.
   max 10*x1 + 15*x2 + 25*x3
   u.d.N. x1 + x2 + 2*x3 <= 40
          2*x1 + 3*x2 + x3 <= 50
          x >= 0

Deckt scipy.optimize, highspy, CVXPY und OR-Tools/GLOP ab, mit Kreuzvergleich
am Ende.

WICHTIG: Jeder Solver laeuft in einem EIGENEN Prozess, weil sich ortools und
highspy auf vielen Systemen nicht gemeinsam importieren lassen (beide bringen
eine eigene HiGHS-Kopie mit -> Symbolkonflikt).

Die Isolation besorgt ein ProcessPoolExecutor. Drei Einstellungen ergeben
zusammen die Garantie:

  mp_context "spawn"     Der Kindprozess startet mit einem FRISCHEN
                         Interpreter, statt den Speicher des Elternprozesses
                         zu erben. Was hier schon importiert ist, ist dort
                         nicht importiert. Mit dem Standard "fork" auf Linux
                         waere das nicht so.
  max_tasks_per_child=1  Jede Aufgabe bekommt einen NEUEN Prozess. Ohne das
                         wuerde der Pool seinen Arbeiter wiederverwenden - und
                         beim zweiten Solver waere der Konflikt zurueck.
  max_workers=1          Haelt die vier Laeufe nacheinander. Nicht aus
                         Vorsicht, sondern damit die gemessenen Zeiten
                         vergleichbar bleiben.

Jeder Solver steht in einer eigenen Funktion mit LOKALEM Import. Das ist der
Unterschied zu einem Codestring, den man an 'python -c' uebergibt: Die
Funktion laesst sich einzeln aufrufen, testen und vom Editor pruefen - ein
String nicht.

Benoetigt: scipy, highspy, cvxpy, ortools
"""

import multiprocessing
import time
from concurrent.futures import ProcessPoolExecutor

ERWARTET = 530.0            # Ergebnis der Handrechnung zum Produktionsprogramm

# Die Instanz - einmal notiert, von allen vier Funktionen benutzt.
ZIEL = [10.0, 15.0, 25.0]
MATRIX = [[1, 1, 2], [2, 3, 1]]
KAPAZITAET = [40.0, 50.0]


def loese_mit_scipy() -> tuple[float, list[float]]:
    from scipy.optimize import linprog
    ergebnis = linprog(c=[-w for w in ZIEL],          # linprog MINIMIERT -> negieren
                       A_ub=MATRIX, b_ub=KAPAZITAET,
                       bounds=[(0, None)] * 3, method="highs")
    return -ergebnis.fun, list(ergebnis.x)


def loese_mit_highspy() -> tuple[float, list[float]]:
    import highspy
    import numpy as np
    h = highspy.Highs()
    h.setOptionValue("output_flag", False)
    h.addVars(3, np.zeros(3), np.full(3, highspy.kHighsInf))
    h.changeObjectiveSense(highspy.ObjSense.kMaximize)
    for j, wert in enumerate(ZIEL):
        h.changeColCost(j, wert)
    # CSR-Format: starts[i] = Beginn von Zeile i in indices/values
    h.addRows(2, np.full(2, -highspy.kHighsInf), np.array(KAPAZITAET), 6,
              np.array([0, 3], dtype=np.int32),
              np.array([0, 1, 2, 0, 1, 2], dtype=np.int32),
              np.array([float(w) for zeile in MATRIX for w in zeile]))
    h.run()
    return (h.getInfo().objective_function_value,
            list(h.getSolution().col_value[:3]))


def loese_mit_cvxpy() -> tuple[float, list[float]]:
    import cvxpy as cp
    import numpy as np
    x = cp.Variable(3, nonneg=True)
    problem = cp.Problem(cp.Maximize(np.array(ZIEL) @ x),
                         [np.array(MATRIX) @ x <= np.array(KAPAZITAET)])
    problem.solve()
    return float(problem.value), [float(v) for v in x.value]


def loese_mit_ortools() -> tuple[float, list[float]]:
    from ortools.linear_solver import pywraplp
    s = pywraplp.Solver.CreateSolver("GLOP")
    x = [s.NumVar(0, s.infinity(), f"x{j+1}") for j in range(3)]
    for i, kapazitaet in enumerate(KAPAZITAET):
        s.Add(sum(MATRIX[i][j] * x[j] for j in range(3)) <= kapazitaet)
    s.Maximize(sum(ZIEL[j] * x[j] for j in range(3)))
    s.Solve()
    return s.Objective().Value(), [v.solution_value() for v in x]


ANSAETZE = {
    "scipy.optimize.linprog": loese_mit_scipy,
    "highspy (natives HiGHS)": loese_mit_highspy,
    "cvxpy": loese_mit_cvxpy,
    "ortools / GLOP": loese_mit_ortools,
}


if __name__ == "__main__":
    print("=" * 78)
    print("  EIN SYSTEM - VIER ANSAETZE (je eigener Prozess)")
    print("=" * 78)
    print(f"{'Bibliothek':<26} {'Z*':>10} {'x1':>7} {'x2':>7} {'x3':>7} {'Zeit':>10}")
    print("-" * 78)

    werte = []
    # Ein Pool, vier Aufgaben, vier frische Prozesse. Der Kontext muss
    # "spawn" sein - siehe Modulkommentar.
    with ProcessPoolExecutor(
            max_workers=1,
            mp_context=multiprocessing.get_context("spawn"),
            max_tasks_per_child=1) as pool:
        for name, funktion in ANSAETZE.items():
            beginn = time.perf_counter()
            try:
                wert, x = pool.submit(funktion).result(timeout=120)
            except Exception as fehler:                     # Bibliothek fehlt o. Ae.
                print(f"{name:<26} nicht verfuegbar: {str(fehler)[:40]}")
                continue
            dauer = time.perf_counter() - beginn
            werte.append(wert)
            print(f"{name:<26} {wert:>10.2f} {x[0]:>7.2f} {x[1]:>7.2f} {x[2]:>7.2f} "
                  f"{dauer:>8.2f} s")

    print("-" * 78)
    spanne = max(werte) - min(werte)
    print(f"Spannweite zwischen den Bibliotheken: {spanne:.2e}")
    print(f"Abweichung zur Handrechnung ({ERWARTET:.0f}):        "
          f"{abs(werte[0] - ERWARTET):.2e}")
    assert spanne < 1e-6, "Die Bibliotheken widersprechen sich!"
    assert abs(werte[0] - ERWARTET) < 1e-6, "Ergebnis weicht von der Handrechnung ab!"
    print("Alle Wege fuehren zum selben, von Hand bestaetigten Optimum.")
    print("(Die Zeiten enthalten Prozessstart und Import - sie messen NICHT die")
    print(" reine Solverleistung. Die Uebungsaufgabe 'Laufzeitvergleich' trennt beides.)")
    print("=" * 78)

modell.kapazitaet = pyo.Constraint(modell.R, rule=kapazitaet)

modell.add_constraints((verbrauch * x).sum("produkt") <= vorrat, name="kapazitaet")

#!/usr/bin/env python3

# Modellierungsschichten.py
"""
Kapitel Oekosystem: Pyomo und Linopy - zwei Modellierungsschichten fuer grosse Modelle.

Geloest wird dasselbe Produktionsproblem wie im Vierfach-Vergleich:

    max 10*x1 + 15*x2 + 25*x3
    u.d.N.  x1 +  x2 + 2*x3 <= 40
           2*x1 + 3*x2 +  x3 <= 50
           x >= 0

Handrechnung: Z* = 530 bei x = (0, 12, 14).

Der Vergleich zeigt die beiden Denkweisen:
  * Pyomo  - algebraisch, indexbasiert, Industriestandard fuer Grossmodelle,
             trennt Modellstruktur sauber von den Daten (AbstractModel).
  * Linopy - beschriftete Arrays (xarray): eine Zeile Code erzeugt Tausende
             Nebenbedingungen auf einmal, ohne Python-Schleife.

Beide bringen KEINEN eigenen Solver mit; hier rechnet in beiden Faellen HiGHS.

WICHTIG: ortools wird in diesem Programm bewusst NICHT importiert - es
vertraegt sich nicht mit der HiGHS-Kopie, die Pyomo und Linopy laden
(siehe die Stolperfalle im Abschnitt 'Ein System - vier Programmieransaetze').

Benoetigt: pyomo, linopy, xarray, pandas, highspy, numpy
"""

from __future__ import annotations

import time

import numpy as np
import pandas as pd
import pyomo.environ as pyo
import xarray as xr
import linopy

ERWARTET = 530.0                     # Ergebnis der Handrechnung

PRODUKTE = ["Standard", "Komfort", "Premium"]
RESSOURCEN = ["Material", "Montage"]

DECKUNGSBEITRAG = np.array([10.0, 15.0, 25.0])
VERBRAUCH = np.array([[1.0, 1.0, 2.0],       # Material je Produkt
                      [2.0, 3.0, 1.0]])      # Montage je Produkt
VORRAT = np.array([40.0, 50.0])


# --- Pyomo: algebraisch und indexbasiert ------------------------------------

def loese_mit_pyomo() -> tuple[float, list[float], float]:
    """Pyomo denkt in Mengen und Indizes, wie ein Mathematiker es aufschreibt.

    `Constraint(RESSOURCEN, rule=...)` erzeugt fuer JEDES Element der Menge
    eine Nebenbedingung - das ist das 'fuer alle i' der Formelsprache,
    unmittelbar in Code uebersetzt.
    """
    t0 = time.perf_counter()

    modell = pyo.ConcreteModel(name="Produktionsprogramm")
    modell.P = pyo.Set(initialize=PRODUKTE)
    modell.R = pyo.Set(initialize=RESSOURCEN)

    modell.db = pyo.Param(modell.P, initialize=dict(zip(PRODUKTE, DECKUNGSBEITRAG)))
    modell.a = pyo.Param(modell.R, modell.P, initialize={
        (r, p): VERBRAUCH[i, j]
        for i, r in enumerate(RESSOURCEN) for j, p in enumerate(PRODUKTE)})
    modell.vorrat = pyo.Param(modell.R, initialize=dict(zip(RESSOURCEN, VORRAT)))

    modell.x = pyo.Var(modell.P, domain=pyo.NonNegativeReals)

    modell.ziel = pyo.Objective(
        expr=sum(modell.db[p] * modell.x[p] for p in modell.P),
        sense=pyo.maximize)

    def kapazitaet(m, r):
        return sum(m.a[r, p] * m.x[p] for p in m.P) <= m.vorrat[r]

    modell.kapazitaet = pyo.Constraint(modell.R, rule=kapazitaet)

    ergebnis = pyo.SolverFactory("appsi_highs").solve(modell)
    dauer = time.perf_counter() - t0

    status = ergebnis.solver.termination_condition
    if status != pyo.TerminationCondition.optimal:
        raise RuntimeError(f"Pyomo meldet Status: {status}")

    return (float(pyo.value(modell.ziel)),
            [float(pyo.value(modell.x[p])) for p in PRODUKTE],
            dauer)


# --- Linopy: beschriftete Arrays --------------------------------------------

def loese_mit_linopy() -> tuple[float, list[float], float]:
    """Linopy denkt in beschrifteten Arrays (xarray).

    Der entscheidende Unterschied: `(verbrauch * x).sum("produkt") <= vorrat`
    ist EINE Zeile und erzeugt so viele Nebenbedingungen, wie die Dimension
    'ressource' Eintraege hat. Bei 2 Ressourcen faellt das nicht auf, bei
    200 000 schon - dort entstehen sie als Matrixoperation statt in einer
    Python-Schleife.
    """
    t0 = time.perf_counter()

    # Benannte Indizes statt blosser Listen: Dadurch heissen die Achsen
    # 'produkt' und 'ressource', und xarray fuehrt sie beim Rechnen von allein
    # richtig zusammen. Ohne Namen vergibt linopy 'dim_0', und man muss
    # spaeter umbenennen - eine haeufige Stolperstelle.
    produkt = pd.Index(PRODUKTE, name="produkt")
    ressource = pd.Index(RESSOURCEN, name="ressource")

    modell = linopy.Model()
    modell.add_variables(lower=0, coords=[produkt], name="menge")
    x = modell.variables["menge"]

    db = xr.DataArray(DECKUNGSBEITRAG, coords=[produkt])
    verbrauch = xr.DataArray(VERBRAUCH, coords=[ressource, produkt])
    vorrat = xr.DataArray(VORRAT, coords=[ressource])

    # EINE Zeile - sie erzeugt so viele Nebenbedingungen, wie die Achse
    # 'ressource' Eintraege hat. Genau das ist der Punkt.
    modell.add_constraints((verbrauch * x).sum("produkt") <= vorrat,
                           name="kapazitaet")
    modell.add_objective((db * x).sum(), sense="max")

    modell.solve(solver_name="highs", output_flag=False)
    dauer = time.perf_counter() - t0

    if modell.termination_condition != "optimal":
        raise RuntimeError(f"Linopy meldet Status: {modell.termination_condition}")

    loesung = modell.variables["menge"].solution.to_series()
    return (float(modell.objective.value),
            [float(loesung[p]) for p in PRODUKTE],
            dauer)


if __name__ == "__main__":
    print("=" * 80)
    print("  MODELLIERUNGSSCHICHTEN FUER GROSSE MODELLE")
    print("=" * 80)
    print(f"{'Schicht':<14} {'Z*':>10} {'Standard':>10} {'Komfort':>10} "
          f"{'Premium':>10} {'Zeit':>10}")
    print("-" * 80)

    ergebnisse = []
    for name, loeser in [("Pyomo", loese_mit_pyomo), ("Linopy", loese_mit_linopy)]:
        ziel, mengen, dauer = loeser()
        ergebnisse.append(ziel)
        print(f"{name:<14} {ziel:>10.2f} {mengen[0]:>10.2f} {mengen[1]:>10.2f} "
              f"{mengen[2]:>10.2f} {dauer:>8.2f} s")

    print("-" * 80)
    for ziel in ergebnisse:
        assert abs(ziel - ERWARTET) < 1e-6, \
            f"Abweichung von der Handrechnung: {ziel} statt {ERWARTET}"
    print(f"Beide stimmen mit der Handrechnung ueberein (Z* = {ERWARTET:.0f}).")
    print()
    print("Wann welche Schicht?")
    print("  Pyomo  -> wenn Modellstruktur und Daten getrennt bleiben sollen,")
    print("            wenn nichtlineare Terme oder MINLP dazukommen koennen,")
    print("            wenn spaeter ein kommerzieller Solver angebunden wird.")
    print("  Linopy -> wenn die Daten ohnehin als beschriftete Arrays vorliegen")
    print("            (Energiesystem-, Netz- und Zeitreihenmodelle) und das")
    print("            Modell zehntausende gleichartige Nebenbedingungen hat.")
    print("=" * 80)

#!/usr/bin/env python3

# Vektorisierte_Modellgenerierung.py
"""
Kapitel Oekosystem: Warum der Solver oft gar nicht der Engpass ist.

In realen Projekten geht ein grosser Teil der Rechenzeit nicht ins Loesen,
sondern ins AUFBAUEN des Modells. Dieses Programm misst das an einem
Transportproblem wachsender Groesse in vier Stufen:

    A  Modellierungsschicht, Nebenbedingung fuer Nebenbedingung (OR-Tools)
    B  Matrix direkt, aber mit Python-Schleifen ueber die Eintraege (COO)
    C  Matrix vektorisiert ueber Kronecker-Produkte (NumPy/SciPy)
    D  Daten kommen als lange Tabelle, aufbereitet mit Polars

Alle Varianten loesen dasselbe Problem und muessen denselben Zielwert
liefern - das wird am Ende geprueft.

Benoetigt: numpy, scipy, ortools; Variante D zusaetzlich polars (optional).
"""

from __future__ import annotations

import importlib.util
import time

import numpy as np
import scipy.sparse as sp
from ortools.linear_solver import pywraplp
from scipy.optimize import linprog

HAT_POLARS = importlib.util.find_spec("polars") is not None


def erzeuge_daten(m: int, n: int, saat: int = 3
                  ) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Transportproblem: m Werke, n Kunden.

    Liefert (kosten[m, n], angebot[m], bedarf[n]). Das Gesamtangebot liegt
    20 % ueber dem Gesamtbedarf, damit das Modell sicher loesbar ist.
    """
    rng = np.random.default_rng(saat)
    kosten = rng.uniform(1.0, 20.0, size=(m, n))
    bedarf = rng.uniform(10.0, 50.0, size=n)
    angebot = np.full(m, 1.2 * bedarf.sum() / m)
    return kosten, angebot, bedarf


# --- Variante A: Modellierungsschicht, Bedingung fuer Bedingung -------------

def loese_mit_modellierungsschicht(kosten: np.ndarray, angebot: np.ndarray,
                                   bedarf: np.ndarray) -> tuple[float, float, float]:
    """So schreibt man ein Transportproblem zuerst hin - gut lesbar, nah an
    der mathematischen Formulierung, jede Nebenbedingung ein eigener Aufruf.

    Jedes `s.Add(sum(...))` baut in Python einen Ausdrucksbaum aus m bzw. n
    Termen auf und uebergibt ihn einzeln an die C++-Schicht. Das ist der
    Preis der Bequemlichkeit - und er waechst linear mit der Modellgroesse.
    """
    m, n = kosten.shape

    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][j] for j in range(n)) <= angebot[i])
    for j in range(n):
        s.Add(sum(x[i][j] for i in range(m)) >= bedarf[j])
    s.Minimize(sum(kosten[i][j] * x[i][j] for i in range(m) for j in range(n)))
    t_aufbau = time.perf_counter() - t0

    t0 = time.perf_counter()
    status = s.Solve()
    t_loesen = time.perf_counter() - t0
    if status != pywraplp.Solver.OPTIMAL:
        raise RuntimeError(f"Solver-Status: {status}")
    return t_aufbau, t_loesen, s.Objective().Value()


# --- Varianten B bis D: Matrix selbst bauen, dann SciPy/HiGHS ---------------

def baue_mit_schleifen(kosten: np.ndarray, angebot: np.ndarray,
                       bedarf: np.ndarray):
    """Die Nebenbedingungsmatrix als COO-Tripel (Zeile, Spalte, Wert), erzeugt
    in verschachtelten Python-Schleifen.

    Schon deutlich naeher am Blech als Variante A - es entsteht kein
    Ausdrucksbaum mehr. Die Schleife selbst bleibt aber Python.
    """
    m, n = kosten.shape
    zeilen: list[int] = []
    spalten: list[int] = []
    werte: list[float] = []
    rechte_seite: list[float] = []

    for i in range(m):                       # Angebot je Werk
        for j in range(n):
            zeilen.append(i)
            spalten.append(i * n + j)
            werte.append(1.0)
        rechte_seite.append(float(angebot[i]))

    for j in range(n):                       # Bedarf je Kunde, als -x <= -bedarf
        for i in range(m):
            zeilen.append(m + j)
            spalten.append(i * n + j)
            werte.append(-1.0)
        rechte_seite.append(-float(bedarf[j]))

    A_ub = sp.csr_matrix((werte, (zeilen, spalten)), shape=(m + n, m * n))
    return A_ub, np.array(rechte_seite), kosten.ravel()


def baue_vektorisiert(kosten: np.ndarray, angebot: np.ndarray,
                      bedarf: np.ndarray):
    """Dieselbe Matrix ohne eine einzige Schleife - ueber Kronecker-Produkte.

    Die Angebotsmatrix ist  kron(I_m, 1_n^T):  je Werk eine Zeile mit Einsen
    an genau den n Spalten dieses Werks. Die Bedarfsmatrix ist
    kron(1_m^T, I_n). Beide entstehen in je einem Aufruf und sind sofort
    duennbesetzt.
    """
    m, n = kosten.shape
    angebots_matrix = sp.kron(sp.identity(m, format="csr"), np.ones((1, n)))
    bedarfs_matrix = sp.kron(np.ones((1, m)), sp.identity(n, format="csr"))

    A_ub = sp.vstack([angebots_matrix, -bedarfs_matrix], format="csr")
    b_ub = np.concatenate([angebot, -bedarf])
    return A_ub, b_ub, kosten.ravel()


def baue_mit_polars(kosten: np.ndarray, angebot: np.ndarray, bedarf: np.ndarray):
    """Der realistische Fall: Die Kosten kommen als LANGE Tabelle
    (werk, kunde, kosten) aus Datenbank, Data Lake oder CSV-Datei.

    Polars berechnet den Spaltenindex jeder Variablen in einem einzigen
    Spaltenausdruck - ohne Python-Schleife ueber die Zeilen. Genau so baut man
    Modelle aus Millionen Tabellenzeilen.
    """
    import polars as pl

    m, n = kosten.shape
    tabelle = pl.DataFrame({
        "werk": np.repeat(np.arange(m), n),
        "kunde": np.tile(np.arange(n), m),
        "kosten": kosten.ravel(),
    }).with_columns(
        (pl.col("werk") * n + pl.col("kunde")).alias("var_index")
    )

    var_index = tabelle["var_index"].to_numpy()
    werk = tabelle["werk"].to_numpy()
    kunde = tabelle["kunde"].to_numpy()
    eins = np.ones(var_index.size)

    angebots_matrix = sp.csr_matrix((eins, (werk, var_index)), shape=(m, m * n))
    bedarfs_matrix = sp.csr_matrix((eins, (kunde, var_index)), shape=(n, m * n))

    A_ub = sp.vstack([angebots_matrix, -bedarfs_matrix], format="csr")
    b_ub = np.concatenate([angebot, -bedarf])
    return A_ub, b_ub, tabelle["kosten"].to_numpy()


def messe_matrixvariante(bauer, kosten, angebot, bedarf
                         ) -> tuple[float, float, float]:
    """Liefert (Aufbauzeit, Loesezeit, Zielwert) fuer die Varianten B bis D."""
    t0 = time.perf_counter()
    A_ub, b_ub, c = bauer(kosten, angebot, bedarf)
    t_aufbau = time.perf_counter() - t0

    t0 = time.perf_counter()
    ergebnis = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=(0, None), method="highs")
    t_loesen = time.perf_counter() - t0

    if not ergebnis.success:
        raise RuntimeError(f"Solver-Status: {ergebnis.message}")
    return t_aufbau, t_loesen, float(ergebnis.fun)


if __name__ == "__main__":
    print("=" * 88)
    print("  MODELLAUFBAU: WO DIE ZEIT WIRKLICH HINGEHT")
    print("=" * 88)
    print("Transportproblem mit m Werken und n Kunden -> m*n Variablen.\n")

    varianten: list[tuple[str, object]] = [
        ("A: OR-Tools, Add() je NB", None),        # Sonderfall, eigener Messpfad
        ("B: COO in Schleifen", baue_mit_schleifen),
        ("C: NumPy vektorisiert", baue_vektorisiert),
    ]
    if HAT_POLARS:
        varianten.append(("D: Polars-Tabelle", baue_mit_polars))
    else:
        print("Hinweis: polars nicht installiert - Variante D wird uebersprungen.\n")

    # Aufwaermlauf: Der erste Aufruf bezahlt Importe und einmalige
    # Initialisierungen. Wer den mitmisst, vergleicht Startkosten statt
    # Rechenarbeit - ein klassischer Benchmark-Fehler.
    aufwaerm = erzeuge_daten(10, 10)
    loese_mit_modellierungsschicht(*aufwaerm)
    for _, bauer in varianten[1:]:
        bauer(*aufwaerm)

    print(f"{'Groesse':<26} {'Variante':<26} {'Aufbau':>9} {'Loesen':>9} "
          f"{'Aufbauanteil':>13}")
    print("-" * 88)

    for m, n in [(40, 40), (120, 120), (250, 250)]:
        kosten, angebot, bedarf = erzeuge_daten(m, n)
        zielwerte = []
        for nummer, (name, bauer) in enumerate(varianten):
            if bauer is None:
                t_aufbau, t_loesen, ziel = loese_mit_modellierungsschicht(
                    kosten, angebot, bedarf)
            else:
                t_aufbau, t_loesen, ziel = messe_matrixvariante(
                    bauer, kosten, angebot, bedarf)
            zielwerte.append(ziel)
            anteil = 100.0 * t_aufbau / (t_aufbau + t_loesen)
            groesse = f"{m}x{n} = {m*n:,} Variablen" if nummer == 0 else ""
            print(f"{groesse:<26} {name:<26} {t_aufbau:>8.3f}s {t_loesen:>8.3f}s "
                  f"{anteil:>12.0f} %")

        # Alle Varianten muessen dasselbe Problem beschreiben.
        spanne = max(zielwerte) - min(zielwerte)
        assert spanne < 1e-6 * max(abs(z) for z in zielwerte), \
            f"Varianten widersprechen sich: {zielwerte}"
        print(f"{'':<26} {'-> Zielwert (alle gleich)':<26} {zielwerte[0]:>9.2f}"
              f"   Spanne {spanne:.1e}")
        print("-" * 88)

    print("\nZwei Lehren aus der Tabelle:")
    print("1. Bei Variante A geht mehr Zeit in den AUFBAU als ins Loesen. Wer hier")
    print("   einen schnelleren Solver kauft, beschleunigt den kleineren Teil.")
    print("2. Zwischen B und C liegt keine andere Mathematik, nur eine andere")
    print("   Schreibweise derselben Matrix - Schleife gegen Kronecker-Produkt.")
    print("\nDie Lesbarkeit von Variante A ist trotzdem viel wert: Fangen Sie dort an,")
    print("und vektorisieren Sie erst, wenn die Messung es verlangt.")
    print("=" * 88)