#!/usr/bin/env python3

# Toleranzen_und_Entartung.py
"""
Kapitel LP: Zwei Faelle, in denen man Schattenpreisen NICHT trauen darf.

  1. Entartung (Degeneriertheit): Mehr Nebenbedingungen sind aktiv, als das
     Problem Variablen hat. Dann ist der Schattenpreis nicht eindeutig - zwei
     korrekte Solver liefern voellig verschiedene Werte fuer dasselbe Optimum.
  2. Toleranzen: "ausgelastet" heisst nie 'schlupf == 0', sondern immer
     'schlupf < toleranz'. Wer auf exakte Gleichheit prueft, baut Berichte,
     die zufaellig mal stimmen und mal nicht.

Teil 3 zeigt die professionelle Antwort auf Fall 1: Statt EINEN Schattenpreis
zu melden, berechnet man seine SPANNE ueber alle optimalen Dualloesungen.

Benoetigt: numpy, scipy
"""

from __future__ import annotations

import numpy as np
from scipy.optimize import linprog

# Ein bewusst entartetes Beispiel: drei Geraden schneiden sich in EINEM Punkt.
#   max x1 + x2
#   u.d.N.  x1 + 2*x2 <= 4        (A)
#          2*x1 +  x2 <= 4        (B)
#           x1 +  x2 <= 8/3       (C)  - laeuft genau durch die Ecke (4/3, 4/3)
# In zwei Dimensionen legen schon zwei Geraden eine Ecke fest. Hier sind drei
# aktiv - eine zu viel. Genau das ist Entartung.
C_ZIEL = np.array([-1.0, -1.0])          # linprog minimiert -> negiert
A_UB = np.array([[1.0, 2.0],
                 [2.0, 1.0],
                 [1.0, 1.0]])
B_UB = np.array([4.0, 4.0, 8.0 / 3.0])
NAMEN = ["A: Fraeszeit", "B: Schleifzeit", "C: Pruefzeit"]


def zeige_entartung() -> float:
    """Loest dasselbe LP mit zwei Verfahren und vergleicht die Dualwerte."""
    print("=" * 78)
    print("  1. ENTARTUNG: DERSELBE PLAN, GEGENSAETZLICHE SCHATTENPREISE")
    print("=" * 78)

    print(f"{'Verfahren':<26} {'x1':>7} {'x2':>7} {'Z*':>9}   Schattenpreise")
    print("-" * 78)

    dualwerte = {}
    for verfahren, beschreibung in [("highs-ds", "Dual Simplex"),
                                    ("highs-ipm", "Innere-Punkte-Verfahren")]:
        ergebnis = linprog(C_ZIEL, A_ub=A_UB, b_ub=B_UB, bounds=(0, None),
                           method=verfahren)
        if not ergebnis.success:
            raise RuntimeError(f"{verfahren}: {ergebnis.message}")
        y = -ergebnis.ineqlin.marginals
        dualwerte[verfahren] = y
        print(f"{beschreibung:<26} {ergebnis.x[0]:>7.3f} {ergebnis.x[1]:>7.3f} "
              f"{-ergebnis.fun:>9.4f}   {np.round(y, 4)}")

    print("-" * 78)
    print("Beide Zeilen sind RICHTIG: gleicher Plan, gleicher Zielwert, und beide")
    print("Dualvektoren erfuellen die Optimalitaetsbedingungen. Trotzdem sagen sie")
    print("das Gegenteil:")
    print(f"  Dual Simplex : {NAMEN[2]} ist wertlos, A und B sind je 0,33 EUR wert.")
    print(f"  Innere Punkte: {NAMEN[0]} und {NAMEN[1]} sind wertlos, C ist 1,00 EUR wert.")
    print()
    print("Wer auf dieser Grundlage eine Maschine kauft, hat eine 50:50-Chance -")
    print("abhaengig davon, welches Verfahren der Solver zufaellig gewaehlt hat.")

    return float(-linprog(C_ZIEL, A_ub=A_UB, b_ub=B_UB, bounds=(0, None)).fun)


def entartung_erkennen() -> None:
    """Der Test, der in jedes Auswertungsskript gehoert."""
    print("\n" + "=" * 78)
    print("  2. ENTARTUNG ERKENNEN - UND WARUM '== 0' DABEI VERSAGT")
    print("=" * 78)

    ergebnis = linprog(C_ZIEL, A_ub=A_UB, b_ub=B_UB, bounds=(0, None))
    schlupf = B_UB - A_UB @ ergebnis.x

    print(f"{'Nebenbedingung':<18} {'Schlupf':>16} {'== 0 ?':>9} "
          f"{'< 1e-7 ?':>10}")
    print("-" * 78)
    for name, s in zip(NAMEN, schlupf):
        print(f"{name:<18} {s:>16.3e} {str(s == 0.0):>9} {str(abs(s) < 1e-7):>10}")

    aktiv = int((np.abs(schlupf) < 1e-7).sum())
    variablen = A_UB.shape[1]
    print("-" * 78)
    print(f"Aktive Nebenbedingungen: {aktiv}, Variablen: {variablen}")
    if aktiv > variablen:
        print(f"=> ENTARTET. {aktiv} aktive Restriktionen bei nur {variablen} "
              "Variablen bedeuten:")
        print("   Der Schattenpreis ist nicht eindeutig. Melden Sie eine Spanne,")
        print("   keinen Einzelwert (siehe Teil 3).")
    else:
        print("=> nicht entartet, die Dualwerte sind eindeutig.")

    print("\nBeachten Sie die Spalte '== 0': Ein Schlupf von 4.44e-16 ist")
    print("rechnerisch null, aber nicht gleich 0.0. Wer mit '==' prueft,")
    print("uebersieht genau die Engpaesse, die er sucht.")


def schattenpreis_spanne(zielwert: float, toleranz: float = 1e-9
                         ) -> list[tuple[float, float]]:
    """Berechnet fuer jede Nebenbedingung die Spanne ihres Schattenpreises
    ueber ALLE optimalen Dualloesungen.

    Die Menge der optimalen Dualloesungen ist selbst ein Polyeder:

        A^T y >= c,   y >= 0,   b^T y = Z*

    (Dualzulaessigkeit plus starker Dualitaetssatz.) Minimiert und maximiert
    man darauf y_i, erhaelt man die exakten Grenzen. Das ist die ehrliche
    Auskunft an das Management: nicht 'die Stunde ist 0,33 EUR wert', sondern
    'zwischen 0,00 und 0,33 EUR - der Wert ist aus dem Modell nicht bestimmbar'.
    """
    m = A_UB.shape[0]
    # A^T y >= c  <=>  -A^T y <= -c ; Ziel war max c^T x, in linprog-Notation
    # steckt c mit negativem Vorzeichen in C_ZIEL.
    c_original = -C_ZIEL
    A_dual_ub = -A_UB.T
    b_dual_ub = -c_original

    spannen = []
    for i in range(m):
        richtung = np.zeros(m)
        richtung[i] = 1.0
        grenzen = []
        for vorzeichen in (1.0, -1.0):          # 1 = minimieren, -1 = maximieren
            ergebnis = linprog(
                vorzeichen * richtung,
                A_ub=A_dual_ub, b_ub=b_dual_ub,
                A_eq=B_UB.reshape(1, -1), b_eq=[zielwert],
                bounds=(0, None))
            if not ergebnis.success:
                raise RuntimeError(f"Spannenberechnung fehlgeschlagen: "
                                   f"{ergebnis.message}")
            grenzen.append(float(ergebnis.x[i]))
        spannen.append((min(grenzen), max(grenzen)))
    return spannen


def zeige_spanne(zielwert: float) -> None:
    print("\n" + "=" * 78)
    print("  3. DIE EHRLICHE AUSKUNFT: SCHATTENPREIS-SPANNEN")
    print("=" * 78)

    spannen = schattenpreis_spanne(zielwert)
    print(f"{'Nebenbedingung':<18} {'von':>10} {'bis':>10}   Aussage")
    print("-" * 78)
    for name, (unten, oben) in zip(NAMEN, spannen):
        if oben - unten < 1e-7:
            aussage = f"eindeutig {oben:.2f} EUR"
        elif oben < 1e-7:
            aussage = "sicher wertlos (kein Engpass)"
        else:
            aussage = "NICHT bestimmbar - Spanne melden!"
        print(f"{name:<18} {unten:>10.4f} {oben:>10.4f}   {aussage}")

    print("-" * 78)
    print("So berichtet man an Entscheider: 'Eine zusaetzliche Fraesstunde ist")
    print("zwischen 0,00 und 0,33 EUR wert - das Modell kann es nicht genauer")
    print("sagen, weil drei Engpaesse exakt gleichzeitig binden.' Das ist eine")
    print("brauchbare Aussage. Ein erfundener Einzelwert ist es nicht.")


if __name__ == "__main__":
    zielwert = zeige_entartung()
    entartung_erkennen()
    zeige_spanne(zielwert)

    print("\n" + "=" * 78)
    print("Merksatz: Pruefen Sie VOR jeder Sensitivitaetsaussage auf Entartung -")
    print("und vergleichen Sie Schlupfwerte nie mit '== 0', sondern mit einer")
    print("Toleranz.")
    print("=" * 78)