operations_research/OR_HTML_04/programme/Benchmark_Skalierung.py

#!/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 json
import subprocess
import sys
import textwrap

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


# Jeder Eintrag ist ein eigenstaendiges Programm: Instanz aufbauen, loesen,
# Ergebnis als JSON ausgeben. Die Instanz wird in jedem Kindprozess aus
# derselben Saat neu erzeugt - so reist nichts ueber die Prozessgrenze,
# was das Ergebnis verfaelschen koennte.
VORSPANN = """
import json, time, resource
import numpy as np

def instanz(m, n):
    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():
    # ru_maxrss ist unter Linux in Kilobyte
    return resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / 1024

M, N = {m}, {n}
kosten, angebot, bedarf = instanz(M, N)
"""

ANSAETZE = {
    "scipy.linprog": """
        from scipy.optimize import linprog
        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
        ausgabe = (float(r.fun), aufbau, loesen, speicher_mb())
    """,

    "highspy": """
        import highspy
        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
        ausgabe = (h.getInfo().objective_function_value, aufbau, loesen, speicher_mb())
    """,

    "ortools/GLOP": """
        from ortools.linear_solver import pywraplp
        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
        ausgabe = (s.Objective().Value(), aufbau, loesen, speicher_mb())
    """,

    "cvxpy": """
        import cvxpy as cp
        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
        ausgabe = (float(problem.value), aufbau, loesen, speicher_mb())
    """,
}


def messe(name: str, quelltext: str, m: int, n: int):
    """Fuehrt einen Ansatz in einem eigenen Prozess aus."""
    programm = (VORSPANN.format(m=m, n=n) + textwrap.dedent(quelltext)
                + "\nprint(json.dumps(ausgabe))\n")
    ergebnis = subprocess.run([sys.executable, "-c", programm],
                              capture_output=True, text=True, timeout=600)
    if ergebnis.returncode != 0:
        return None, ergebnis.stderr.strip().splitlines()[-1][:60]
    return json.loads(ergebnis.stdout.strip().splitlines()[-1]), None


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, quelltext in ANSAETZE.items():
            werte, fehler = messe(name, quelltext, 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)